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ABSTRACT 


Mathematical models of aircraft dynamics typically contain quantities 
called parameters, which depend, in general, on the flight condition and the 
aircraft geometry. It is important to be able to estimate these parameters 
accurately from flight testing the aircraft. The parameter estimates from flight 
data are used to corroborate wind tunnel parameter estimates, validate and 
improve a priori aerodynamic calculations, update aircraft dynamic models 
for flight control system refinement, and predict aircraft responses for realistic 
flight simulation and ground-based pilot training. 

Today, many fighter aircraft are designed with inherent longitudinal 
static instability in order to enhance maneuverability and performance. 

These aircraft must employ stability augmentation systems which, because of 
safety considerations, cannot be turned off. In addition, multiple control 
surfaces are required for enhanced performance in expanded flight envelopes, 
such as high angle of attack flight. 

The object of this research is to develop an algorithm for the design of 
practical, optimal flight test inputs for aircraft parameter estimation 
experiments. This algorithm must be capable of designing multiple input 
experiments for estimation of open loop model parameters from (necessarily) 
closed loop flight test data. 

A general, single pass technique was developed which allows global 
optimization of the flight test input design for parameter estimation using 
the principles of dynamic programming. Provision was made for practical 
constraints on the input form, including amplitude constraints, control 
system dynamics, and selected input frequency range exclusions. In addition, 
the input design was accomplished while imposing output amplitude 
constraints required by model validity and considerations of safety during the 
flight test. The algorithm has multiple input design capability, with optional 
inclusion of a constraint that only one control move at a time, so that a 
human pilot can implement the inputs. The dissertation includes a new 
formulation of the optimal input design problem, a description of a new 
approach to the solution, and a summary of the characteristics of the 
algorithm, followed by three example applications of the new technique 
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which demonstrate the quality and expanded capabilities of the input designs 
produced by the new technique. In all cases, the new input design approach 
showed significant improvement over previous input design methods in 
terms of achievable parameter accuracies. 

The work described in this document was done in partial fulfillment of 
the requirements for the Doctor of Science degree in Aerospace Engineering at 
The George Washington University Joint Institute for the Advancement of 
Flight Sciences (JIAFS), NASA Langley Research Center, Hampton, Virginia. 
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Chapter I - Introduction 


Aircraft flight tests designed specifically for parameter estimation are 
generally motivated by one or more of the following objectives: 

1. The desire to correlate aircraft aerodynamic characteristics obtained 
from wind tunnel experiments and aerodynamic calculations with flight 
test data. 

2. Refinement of the aircraft model for control system analysis and design. 

3. Accurate prediction of the aircraft response using the mathematical 
model, including flight simulation and flight envelope expansion. 

4. Aircraft acceptance testing. 

The design of an experiment to achieve any of the above objectives 
involves specification of the instrumentation and signal conditioning, the 
flight test operational procedure, the inputs for the flight test maneuver, the 
model structure, and the parameter estimation algorithm. In this 
dissertation, the maneuver design, or, specifically, the design of flight test 
input signals, will be studied independently of the other aspects which impact 
the success of the flight test , Other considerations in the flight test design can 
be accounted for in the detail of the input design problem formulation. The 
fundamental principles and procedures regarding the design of the input 
remain unaltered. 

In order to obtain the most accurate estimates of aircraft model 
parameters, the information content in the aircraft response during the flight 
test must be maximized. In general, an aircraft model contains multiple 
response variables and multiple aircraft model parameters. The information 
contained in the aircraft response is embodied in a matrix called the 
information matrix, whose elements are combinations of partial derivatives 
of the aircraft response variables with respect to the model parameters. These 
partial derivatives are called sensitivities, and are obtained by solving the 
so-called sensitivity equations that result from differentiating the state and 
output equations for the aircraft model with respect to each of the parameters 
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in the model. Information matrix elements also depend on the 
measurement sampling rate, and the measurement noise characteristics, 
which are indirectly specified when selecting the instrumentation system. 

In this work, the model structure is assumed known, and aircraft 
model parameters are assumed to be estimated from the flight test data using 
an asymptotically unbiased and efficient parameter estimation technique 
known as maximum likelihood estimation [1]. It is not necessary that the 
parameter estimation algorithm be specifically maximum likelihood. Any 
other asymptotically unbiased and efficient estimator could be used instead. 
The accuracies of the estimated model parameters are given by parameter 
standard errors, which are computed as part of the parameter estimation 
algorithm. The standard error is the value of one standard deviation 
associated with the estimate of a model parameter. It can be shown that the 
theoretical best (lowest) values for the parameter standard errors depend only 
on the information content of the experiment, as embodied in the elements 
of the information matrix. These theoretical lower bounds on the parameter 
standard errors are referred to as Cramer-Rao bounds. The Cramer-Rao 
bounds are independent of the parameter estimation algorithm used to 
extract parameter estimates and standard errors from the data records, 
provided that the parameter estimation algorithm is asymptotically unbiased 
and efficient. Thus, the merit of an input design for aircraft parameter 
estimation can be determined by examining the Cramer-Rao bounds, since 
the latter depend only on the information matrix resulting from the response 
of the aircraft to the application of the input. Comparisons using the 
Cramer-Rao bounds separate the merits of the input design from the merits 
of the parameter estimation algorithm used to extract the model parameter 
estimates from the data. In other words, input designs are compared based 
only on the information content in the experiment, the latter being calculable 
before any parameter estimation is done, and thus independently of the 
particular algorithm used to estimate the values of the parameters. 

The performance, and thus the optimality, of an input design depends 
on the values of the Cramer-Rao bounds associated with that input. Implicit 
in the computation of the Cramer-Rao bounds is the a priori dynamic system 
and measurement model. When designing inputs for aircraft parameter 
estimation flight tests, there must be a complete model ( including parameter 
values ) of the physical system to use during the process of the input design. 
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That is, a model complete with parameter values is necessary in order to 
design an experiment which will produce estimates of the parameters in the 
model. This has been called the "circularity" problem [2]. For aircraft, the 
problem is mitigated by use of parameter estimates obtained from either 
aerodynamic calculations or wind tunnel experiments. The calculation of the 
Cramer-Rao bounds also implicitly includes the time length of the flight test 
data run, along with any imposed constraints on either the input form or 
output response variables. 

In the past, research works which addressed the input design problem 
for airplane parameter estimation used a linear time-invariant dynamic 
model, where the state variables were perturbations about some trim 
condition [2] [3] [4] [5] [6] [7] [8]. The most prevalent approach to optimizing 
the input for aircraft parameter estimation experiments was first done by 
Mehra [5]. The problem was formulated as a fixed time optimal control 
problem, with minimization of a time integral of a scalar function of the 
information matrix as the criterion of optimality. Some difficulty has arisen 
as a result of the desire to find the best input with respect to a scalar measure 
of optimality. A scalar norm of the information matrix must be used to 
encompass the multiple sensitivities contained in the information matrix. 
Denoting the information matrix by M, previous works have maximized 
Tr[M], or det [M], or minimized Tr[M ] ], det [M _1 ], or X max [M T ] where X max is 
the maximum eigenvalue of M" 1 , among others. There has been some debate 
as to which matrix norm should be used to obtain the "best” set of parameter 
standard errors. Reference [6] gives a discussion and review of the different 
matrix norms used by various researchers. Regardless of the scalar optimality 
criterion used for the input optimization, comparisons regarding the 
performance of the designed inputs must always be made by examining the 
Cramer-Rao bounds, since these are the theoretical lowest parameter standard 
errors achievable with a given input. The validity of using the Cramer-Rao 
bounds to assess relative efficacy of input designs for airplane parameter 
estimation experiments has been verified experimentally [2] [7]. The 
Cramer-Rao bounds are of principal significance in designing inputs for 
parameter estimation experiments. 

In practice, there are often subsets of the entire parameter set which are 
of greater importance, and hence need to be estimated with greater accuracy 
than the remainder of the parameter set. The subsets depend strongly on the 
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end purpose of the parameter estimates. This situation considerably 
complicates the problem of choosing a suitable matrix norm of M for the 
optimal input design. In past work, a constant weighting matrix has been 
introduced in association with the information matrix, M, in order to account 
for varying importance of subsets of the entire parameter set [6]. 

Researchers who used any of the various matrix norms of M to design 
inputs for aircraft parameter estimation also used a fixed test time [2] [3] [4] [5] 
[6] [7] [8], The value of the fixed test time was usually chosen by a heuristic 
argument based on previous experience. Chen [9] was the first to realize that 
the true situation in practice is that specific goals for the Cramer-Rao bounds 
could be specified a priori, and that the appropriate goal of the flight test was 
not to optimize some norm of M over a fixed time, but to reduce the ( limited 
and expensive ) flight test time required to achieve desired goals for the 
Cramer-Rao bounds. This approach eliminates the need for any weighting 
matrices, since the objective has been changed to a direct requirement that 
Cramer-Rao bounds be less than or equal to specified target values. 

Parameters with more stringent accuracy requirements are simply assigned 
lower target values for their respective Cramer-Rao bounds. This procedure 
implicitly implements parameter weighting. A priori Cramer-Rao bound 
goals are a function of the purpose for which the parameter estimate values 
are intended, and can be specified by the end user of the parameter estimates. 
Unfortunately, the Chen approach required extensive iteration and 
considerable judgment on the part of the analyst for its use. Rather than use 
an optimization procedure, Chen simply tried a number of candidate input 
designs, which were generated as members of an orthogonal function set 
called Walsh functions. Thus, the Chen solution lacked optimality 
properties, in that no conditions which correspond to an optimal solution 
(such as a zero gradient of the cost with respect to the input vector) were 
satisfied. 

In many previous works, the input design for aircraft parameter 
estimation was done with the input subject to an energy constraint consisting 
of a time integral of the square of the input amplitude. Values used for the 
energy constraint were presumably based on experience with similar flight 
test situations. In cases where the optimal input problem has been 
formulated as an optimal control problem using a scalar norm of the 
information matrix integrated over a fixed time as the optimality criterion, an 
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energy constraint on the input was preferred due mostly to the computation 
simplicity which results from such a constraint. The justification of this 
constraint form has been that the designed input amplitudes are kept at a 
reasonable level, since the test time has been fixed, and the resulting output 
amplitudes are also (indirectly) limited. Another type of input energy 
constraint was used in the work of Reid [3], who required that the input 
amplitude be full positive or full negative amplitude at any time during a 
fixed test time. As shown in Chapter II, this amounts to a constraint on the 
input energy. 

The relatively small number of references for this dissertation reflects 
the fact that much of the research done on optimal input design addresses 
energy constrained, single input, single output problems, often in the 
frequency domain, and often with few (one or two) model parameters. These 
works were not considered relevant to the general multiple input, multiple 
output, multiple parameter, relatively short time, optimal input design 
problem for aircraft parameter estimation, and were therefore omitted from 
the reference list. To the author's knowledge, the reference list given here 
represents the significant work in optimal input design for aircraft parameter 
estimation experiments. 

The main purpose of this dissertation is to describe and demonstrate a 
new approach to the optimal input design problem for aircraft parameter 
estimation experiments. Motivation for the development of the new 
technique arose from several important considerations in the modern flight 
test environment which were either addressed poorly or not at all in 
previous works. These considerations are briefly outlined below. 

Multiple input design - Flight testing of modern aircraft for parameter 
estimation often requires a multiple input design. Here, "input" can mean 
either control surface deflections directly, or pilot station inputs, depending 
on the flight test equipment. Modern aircraft have numerous control 
surfaces in order to achieve expanded flight envelopes and capabilities, and to 
improve failure robustness. The work described here provides multiple 
input design capability in a straightforward way. 

Practical output constraints - Two important considerations require that 
inputs for parameter estimation flight experiments be designed so that output 
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variable amplitudes remain within specified limits. First, the safety of the 
pilot and the aircraft requires that output variable amplitudes be constrained, 
so as to avoid unusual aircraft attitudes from which the pilot cannot recover, 
or excessive accelerations or rotation rates which may damage the aircraft. 
Second, since the aircraft model is typically valid only over certain ranges of 
amplitudes of certain output variables, an experiment which takes the aircraft 
outside the flight regime where the mathematical model is assumed to be 
valid degrades the quality of the data for parameter estimation purposes. 
Previous works relied on the input energy constraint to keep output variable 
amplitudes reasonable. This approach becomes particularly difficult to use for 
multiple inputs or closed loop models, along with being unrepresentative of 
the true practical flight test situation. The present approach includes output 
amplitude constraints directly as part of the optimal input design problem 
formulation. 

Cosed loop models - Many modern aircraft are designed to be open loop 
unstable in order to enhance performance. This requires the use of full time 
stability augmentation systems (SAS) which employ automatic closed loop 
control. These systems cannot be turned off for any length of time because of 
safety considerations. For this case, it is necessary to include a model of the 
control system when analyzing the aircraft response, in order to be able to 
separate the effects of the automatic control system from the open loop 
aerodynamics of the airframe. The aircraft open loop parameters are of 
principal interest, since these are required to realize the objectives stated at 
the beginning of this chapter. The technique described here can accommodate 
the model structure necessary to design inputs for estimating open loop 
model parameters under the condition that only closed loop response data are 
available. The system response must, however, be limited to a region of 
linear control effectiveness and linear open loop dynamics. When 
nonlinearities are involved, the appropriate model structure is a function of 
the particular input used to excite the system [6]. But, the model structure is 
needed in the first place in order to design the input for parameter estimation 
purposes. Thus, input design for nonlinear systems involves a complex 
circularity. The problem of input design for nonlinear models is recognized 
as important and worthy of attention; however, this problem is also 
considered outside the scope of the present investigation, and is not 
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addressed. It was assumed that the capability for imposing arbitrary output 
amplitude constraints could be used in any flight condition to restrict the 
system behavior to a regime which could be adequately described by a linear 
model structure. 

Practical control implementation - Regardless of whether a human pilot or a 
computer implements the designed inputs, it must be assured that the inputs 
can be realized. This means that control system dynamics should be included 
in the input design analysis insofar as possible. When a human pilot must 
implement the input design, allowance must be made for the limited 
accuracy and repeatability of input forms generated by humans, particularly 
for the case of multiple inputs. Additional practical limitations pertain to the 
input amplitude and the input frequency spectrum. It is preferable to include 
these practical constraints on the input form as an integral part of the analysis 
and input design procedure. The present approach does this. 

Minimum flight test time - Resources are always limited; therefore, the flight 
test experiment should be designed to minimize the flight test time required 
to answer the questions which engendered the idea of an experiment. 
Minimizing the flight test time in the context of the present work is 
equivalent to designing so that the a priori Cramer-Rao bound goals can be 
achieved in the minimum flight test time. This factor is significant due to 
the limited availability of flight test aircraft, and the great expense associated 
with flight testing modern aircraft. In addition, for the case of flight testing at 
high angle of attack, the time available for flight maneuvers at a given flight 
condition is limited by altitude loss during the flight test. The approach 
described here uses the principles of dynamic programming and Bellman's 
principle of optimality [10] [11] to produce an input which is globally time 
optimal, subject to the constraints of the problem formulation. 

Maximum parameter accuracy for a fixed flight test time - In practical flight 
test situations, it is often necessary to design an input for a fixed flight test 
time, chosen a priori. The principles of dynamic programming can be applied 
to this situation as well, rendering a global optimal solution for a chosen cost 
function which incorporates the minimization of the Cramer-Rao bounds, 
again subject to the constraints of the problem formulation. 
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This work is an exposition of a new approach to the optimal input 
design problem for aircraft parameter estimation experiments, as well as a 
demonstration of the capabilities of the new technique for producing 
practical, optimal input designs. The remainder of the dissertation will 
describe the problem formulation and solution methodology, followed by a 
summary of the solution algorithm characteristics and three example 
applications of the new technique. The first example is the solution of a 
problem studied by other researchers, and demonstrates improved input 
designs using the new approach. The second example is the solution of a 
problem not treated elsewhere, which highlights the multiple input 
capability, output amplitude constraint capability, and the practical input 
form constraint features included in the new technique. Finally, the third 
example demonstrates the optimal input design technique as applied to a 
closed loop model of the F-18 fighter aircraft longitudinal dynamics. A six 
degree of freedom nonlinear simulation which uses tabular wind tunnel data 
for the vehicle aerodynamics was used as the test aircraft. A summary and 
conclusions section completes the dissertation. 
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Chapter 11 - Problem Formulation 


The goal of the work described here is to design optimal flight test 
inputs for a parameter estimation experiment. Optimality is a property 
whose meaning depends on the problem formulation. In general, the 
problem is specified by the dynamic system and measurement model, the 
criterion of optimality, the admissible control set, and the applicable 
constraints. In what follows, each of these aspects of the problem formulation 
will be discussed regarding its relationship to the physical situation of flight 
test for parameter estimation. 

The dynamic system and output models are linear, with an 
n-dimensional state vector, x, and system dynamic, control and observation 
matrices F, G, and H, respectively, which, in general, depend on a 
p-dimensional parameter vector, 0. The state vector has zero initial 
conditions. The latter condition indicates the assumption of a perturbation 
model derived relative to some specified condition. This approach provides 
known initial conditions, which may be taken as zero with no loss of 
generality. An input vector, u, of dimension m is assumed. The linear 
output model has an output vector y of dimension q, which depends on x and 
0 . Amplitude constraints are imposed on all inputs and on selected outputs, 
and are denoted by pj and respectively. The measurements y m are made at 
N sampling times, separated by constant sampling intervals of length At, with 
additive white Gaussian measurement noise, u. The statistical properties of 
the measurement noise are assumed known ( e g., from ground calibration of 
the instrumentation system). The preceding is expressed mathematically as: 


x(l) = F (0) x ( t) + G (0) u ( t) , x(0) = 0 

y(t) = H (0) x(t) 

y m (i) = y(i) + U(i) i = 1, 2, 3, ..., N 


( 2 - 1 ) 

( 2 - 2 ) 

(2-3) 
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The measurement noise u(i) is assumed Gaussian with 


E{v(i)}=0 , E {\)(i)u T (j)J =R6 jj (2-4) 

where E{ • } is the expectation operator, and 8jj is the Kronecker delta. 

The following constraints are imposed on input and output 
amplitudes: 


1 Uj (t) | < Hj V t, j=l, 2, ...,m 

(2-5) 

I y k (t> | £ n k V t/ ke (1,2, ...,q) 

(2-6) 


Equations (2-1) through (2-6) give the general form of the dynamic 
system and measurement model used in this work. A linear dynamic and 
aerodynamic model with multiple inputs is assumed, with amplitude 
constraints imposed on all inputs and selected outputs. 

The issue of the robustness of the input design to inaccuracies in the a 
priori values for the model parameters is not addressed in this work. This 
problem is important and complex, but cannot be adequately treated here. 

For the simple case of a single input, single output model with one 
model parameter, the best input for a parameter estimation experiment is the 
input which maximizes the sensitivity of the output quantity to changes in 
the parameter over the test time, T. This can be expressed as 


-.-l 


max 

u* t) = 

N 

I 

f dy(i) j 

min 

N 

I 

dy(i)' 

2~ 

u(t) e U 

i = l 

L ae / 

~ u(t) e U 

i=l 

j 50 j 



(2-7) 


where u*(t) is the scalar optimal control, U is the set of admissible controls, 
and the summation over N time points approximates a time integral, with 

T = N At (2-8) 


where At represents the sample time. 


10 



The optimization of an input with respect to the scalar criterion in 
equation (2-7) is straightforward. In most practical situations of interest, 
however, there are multiple outputs and multiple model parameters. In this 
case, the sensitivity of the outputs to changes in the model parameters is 
embodied in a matrix called the information matrix. As discussed in Chapter 
I, in this case it is not clear what scalar optimization criterion should be used 
for the optimal input design. The (p x p) information matrix, M, is given y 
the following expression, derived in Appendix A: 


M = 



(2-9) 


The (q x p) matrix of sensitivities in the above expression is sometimes 
referred to as the discrete sensitivity matrix. Si : 

S _ 3y(i) 

' 30 ( 2 - 10 ) 


The partial derivatives which appear in the expression for the information 
matrix are found using the so-called output sensitivity equations, which are 
obtained by differentiating (2-2) with respect to the components of the 
parameter vector, 0, denoted by 0 k , 

iy_ = H — + ~x k=l,2,...,p 

30 k 30 k 00 k (2-11) 


The partial derivatives of the state vector, x, with respect to 6 k , known as the 
state sensitivities, are found by differentiating (2-1) and switching the order of 
the differentiation on the left side of the equation ( the latter procedure is 
valid as long as x is analytic, which is assumed ) to obtain 


d 

dx 

dt 

30, 



3F 

30 k 


x + 


3G 

u 

30 k 


k=l,2,...,p 

(2-12) 
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The (n p) state sensitivity equations (2-12) must be solved, along with the (n) 
system dynamic equations (2-1), with the results used in the (q p) output 
sensitivity equations (2-11), in order to assemble the information matrix from 
equation (2-9). All initial conditions for equations (2-12) are zero, since the 
initial state is assumed to be known. Thus, 


3x(0) 

"ae^ 


= o 


k = l,2,...,p 


(2-13) 


The information matrix depends on the input indirectly through the 
output and state sensitivity equations, (2-11) and (2-12), because the state 
sensitivity dynamics depend on the input u directly, and both the output 
sensitivities and the state sensitivity dynamics depend on x, which is a 
function of u by the dynamic equations (2-1). 

In Appendix A, it is shown that the inverse of the information matrix 
is the theoretical lower limit for model parameter covariances computed 
using any asymptotically unbiased and efficient parameter estimation 
algorithm with the flight test data used to assemble the information matrix. 
This theoretical lower limit is often called the dispersion matrix, D, and is 
given by 


D = M 


Using (2-9), equation (2-14) may also be written as 


D = 


i-i\ 39 I \ 39 I 


'W 


or, using (2-10), 


D = 


-.1 


N t 

I Si R Si 

i = 1 


(2-14) 


(2-15) 


(2-16) 
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M and D are symmetric, positive definite matrices, since R must be symmetric 
and positive definite, and there can be no zero rows or columns in Si, as long 
as the model structure includes all the parameters. 

M depends on the input time history; therefore, it follows from (2-14) 
that the theoretical lower bounds on the parameter covariances are also 
functions of the input time history. The diagonal elements of the dispersion 
matrix, D, are the theoretical minimum values of the individual parameter 
variances. The Cramer-Rao lower bounds for the parameter standard errors 
are given by the square root of the diagonal elements of the dispersion matrix. 
These values are denoted by a k , 

a k = ^kk / k=l,2,...,p ( 2-17) 

where dj k are the matrix elements of the dispersion matrix, D, i.e., 

^ “ [ ^jk ] > j = 1/2,.. .,p , k = l,2,...,p (2-18) 

The k th diagonal element of the dispersion matrix corresponds to the k*h 
parameter in equations (2-11) and (2-12). The o k values in equation (2-17) may 
be thought of as the minimum possible value of one standard deviation (the 
standard error) for the estimate of the k*h parameter. 

The present work borrows from Chen's [9] view of the optimal input 
design problem, in that the flight test time is minimized for attaining 
Cramer-Rao bound goals which are specified a priori. From equations (2-15) 
and (2-17), the achievable Cramer-Rao bounds from the dispersion matrix are 
a function of the total test time, T, since T = NAt by equation (2-8). As a result 
of this, an inferior input design will require a longer test time,T, to achieve 
specified Cramer-Rao bound goals. 

There are actually two optimal input design cases solved in this work, 
both of which are useful for practical input designs in realistic flight test 
situations. The two cases are the minimum time optimal input design and 
the fixed time optimal input design. 

In the first case, the optimal input is that which achieves specified goals 
for the Cramer-Rao bounds in minimum time , subject to the constraints of 
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the problem formulation. This is a minimum time problem, so that the cost, 
J, is given by 


J = T when o k Vk=l,2,...,p (2-19) 

where £k are the goal values for the Cramer-Rao bounds specified a priori, 
and T is now a variable flight test time to be minimized. Minimizing J in 
equation (2-19) is equivalent to finding the input time history with 
minimum N in equation (2-15) when the condition (dkk) 1 ^ 2 — Ck V k=l,2,...,p 
is satisfied. 

In the second case, the flight test time, T, is fixed, and the optimal input 
is that which achieves minimum Cramer-Rao bounds over a fixed time, T, 
subject to the constraints of the problem formulation. This is a fixed time 
problem, with the cost function 


j = tK ) 2 

k=l (2-20) 

where <Jk are the Cramer-Rao bounds from equation (2-17), and N is fixed in 
equation (2-15). 

In chapter III, the method of solution is outlined for both the 
minimum time problem, with cost given by equation (2-19), and the fixed 
time problem, with cost given by equation (2-20). The topic of constraints is 

addressed next. 

In past work, optimal input design for aircraft parameter estimation 
was carried out with the input subject to an energy constraint of the form: 

C T 

f u(t) T u(t) dt = E 

Jo (2-21) 

where E is some fixed value of the allowable input energy. The interval of 
integration, [ 0,T ], was fixed also, and corresponds to the flight test time. 
Values used for the energy constraint, E, and the flight test time, T, were 
typically chosen by heuristic arguments based on experience with similar 

flight test situations. 
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In practice, there is no direct constraint on the amount of input energy 
which can be applied during the flight test, since neither the pilot nor any 
control system have inherent energy limitations. The practical flight test 
situation dictates that the constraints be directly on the amplitudes of both the 
input and the output variables, as given by equations (2-5) and (2-6), 
respectively. Input amplitude constraints are necessary because any aircraft 
has limitations on control surface movements due to mechanical stops, flight 
control software limiters, or linear control effectiveness limitations. In 
addition, the safety of the aircraft and pilot during the flight test, as well as the 
validity of the aircraft model, require that selected output variable amplitudes 
stay below threshold values throughout the flight test. The constraint form 
(2-21) attempts to achieve all of these practical amplitude constraints in an 
indirect and imprecise way, in order to reap the benefit of computational 
simplicity for a fixed time optimal control solution using variational calculus. 
Use of the constraint (2-21) compromises the effectiveness of the input by not 
allowing full advantage to be realized within whatever amplitude constraints 
are imposed on the input and output variables by the circumstances of the 
flight test. Previous studies have imposed the constraint (2-21), and then, 
after the design was completed, checked the input and output amplitudes 
through simulation to make sure that the practical amplitude constraints 
were not violated during the flight test [2] [3] [4] [5] [7] [8] [9] [12]. The present 
work solves the optimal input design problem with amplitude constraints 
imposed directly on the input and output variables as part of the problem 
formulation. Further background material and rationale for this decision 
follows. 

Define a "bang-bang" input as an input whose amplitude can only be 
full positive or full negative at any time, which means that the admissible 
control set, Ub, has only two elements for a scalar input, namely full positive 
amplitude and full negative amplitude. 


u b=( +H' -V )' 


( 2 - 22 ) 


u(t)eU b , t e [0,T] 


(2-23) 
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where +|i and -g are the upper and lower input amplitude constraints, 
respectively. 

For any "bang-bang" input form over a fixed time interval, the energy 
constraint (2-21) is equivalent to an amplitude constraint. For example, the 
"bang-bang" input signals shown in figures 1(a) and 1(d) have the same total 
test time, and different energy, as computed from equation (2-21). It is thus 
necessary that these signals also have different amplitudes, as can be seen 
from figures 1(a) and 1(d). In fact, for these fixed time "bang-bang" input 
forms, energy and amplitude constraints are related by 

E = p 2 T (2-24) 

For input forms other than the "bang-bang" type, this equivalence of 
amplitude and energy constraints is not valid; however, in some approximate 
sense, it may be claimed that the constraint in equation (2-21) amounts to an 
amplitude constraint on the input when the test time is fixed. The amplitude 
constraint via an energy constraint is imprecise to the extent that the input 
form being considered differs from the "bang-bang" type of input form. 

In past work, various input designs were compared to one another by 
designing inputs with the same energy, computed from equation (2-22), and 
then comparing the resulting diagonal elements of the dispersion matrix 
resulting from each designed input. As seen from figures 1(a), 1(b), and 1(c), 
"bang-bang" input forms with the same energy have different amplitudes 
when the total test time is different. A review of equations (2-9) through 
(2-18) reveals that the Cramer-Rao bounds from the diagonal of the 
dispersion matrix are nonlinear functions of the input time history. The 
preceding statement is true even in the case of linear dynamic and output 
models, by virtue of equation (2-15). It follows in general that the diagonal 
elements of the dispersion matrix are nonlinear functions of the input 
amplitude. Thus, the overall performance of the input design is a nonlinear 
function of the input amplitude. This fact is highlighted when considering 
the different responses of a dynamic system to a high amplitude, short 
duration input such as figure 1 (b), as opposed to a low amplitude, long 
duration input like figure 1(c), which have the same energy, computed from 
equation (2-21). It is clear that amplitude and test time may be adjusted over 
wide ranges without changing the input energy, while at the same time 
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drastically changing the system responses on which the merit of an input 
design is based. The diagonal elements of the dispersion matrix (ok from 
equation (2-17)) are the "bottom line", so to speak, when comparing various 
input designs for parameter estimation experiments. The above 
considerations raise questions regarding the validity of any previous 
comparisons done among input forms where only the input energy was kept 
constant. This type of comparison has been done routinely in the 
literature [2] [4] [6] [7] [9] [12]. 

Several interrelated attributes characterize an input design. This 
renders the task of equitable comparison among various input designs 
problematic. The above discussion concerning figure 1 pointed up the 
relationship between input amplitude, input energy, and total test time. In 
light of that discussion, the actual energy constraint on the input form is a 
composite result of the practical constraints on the input amplitude, the 
control system dynamics, and the flight test time. Thus, if one assumes the 
viewpoint that the flight test time should be minimized, there is no way to 
specify a priori what the value of this input energy constraint should be. 

Even with the test time and input amplitude fixed, control system dynamics 
can have a significant impact on the value of the input energy, so that it is 
still difficult to specify the value of the input energy constraint E in equation 
(2-21), a priori. 

In general, the input form of figure 1(b) is often more effective in the 
sense of achieving specific Cramer-Rao bound goals in minimum time, 
compared to the input form in figure 1(c). This is the result of a more nearly 
impulsive input to the dynamic system in figure 1(b), even though the input 
energy is identical to that in figure 1(c). Thus there exists an interplay 
between allowable amplitude and required test time for given Cramer-Rao 
bound goals -- higher amplitudes result in shorter test times and lower 
amplitudes result in longer test times. For fixed time problems, higher input 
amplitude constraints are associated with potentially lower Cramer-Rao 
bounds for the same flight test time. 

There are further interrelations. Assume now that a fixed output 
amplitude constraint is imposed, and consider the frequency of the input 
design in figure 1(a). In order to satisfy the fixed output amplitude constraint, 
the allowable input amplitude is related to the frequency of the input signal. 
For example, a low frequency input will require a relatively low amplitude so 


17 



I 


that the resulting system output response will remain inside the fixed output 
amplitude constraint. Higher input frequencies would allow larger input 
amplitudes. But the dynamic system output response magnitudes depend on 
input frequency regardless of any constraints on the output amplitudes, since 
different input frequencies excite different dynamic response modes. This 
again impacts the allowable input amplitude as a function of the frequency of 
the input signal. But these effects on input amplitude also impact the 
required test time for minimum time solutions, or the achievable parameter 
accuracy for fixed time solutions, as mentioned above. The interrelations 
among input characteristics are seen to be complex, even for the simple 
"bang-bang" input forms of figure 1. The complexity arises from the fact that 
the dynamic system and output equations are involved, along with the 
matrix inversion of equation (2-15). Input forms other than "bang-bang" (e.g., 
sinusoidal) are subject to the same arguments in general; however, the issues 
are more difficult to discern. 

In addition to the foregoing, there are practical considerations 
encountered when attempting to implement an input design in a real flight 
test environment. One such consideration is that the allowable frequencies 
contained in the input time history will be limited at the high end, due to 
limited instrument dynamic response, reduced dynamic system response to 
high frequency input, and the high frequency limitations of the pilot and 
control system. Another frequently overlooked issue might be called 
implementation distortion. This refers to the practical fact that a designed 
input may be distorted when actually implemented, owing to the limited 
capabilities of the pilot, or to control system dynamics, such as lag and 
hysteresis. As might be expected, some input designs are less susceptible than 
others to these practical difficulties. Accounting for these difficulties during 
the problem formulation can improve the input design. The optimal input 
design procedure described in the next chapter successfully addresses these 
and other constraint issues. 

Previous studies of optimal input design for aircraft parameter 
estimation treated the topic of output constraints superficially. In some cases, 
the energy constraint on the input given in equation (2-21) was claimed as an 
indirect constraint on the output amplitudes. It has been argued that this type 
of indirect constraint on the output is appropriate, due to the gradual 
degradation of the descriptive capability of linear models with increasing 
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excursions from the trim condition [2]. However, this gradual degradation 
cannot be relied on for arbitrary flight conditions. Even if a gradual 
degradation could be guaranteed, such an approach cannot precisely limit 
output amplitudes, and thus compromises the effectiveness of the input by 
possibly not taking full advantage of the system response potential, subject to 
the practical limits on the output amplitudes. Typically in past studies, a 
design was completed, and the aircraft response was checked after the fact 
through simulation to assure that the output amplitudes did not exceed the 
threshold values. The direct method of including output amplitude 
constraints in the problem formulation has been alluded to in the literature, 
but not yet demonstrated, apparently because of the added computational 
complexity which would be introduced. 

For a real flight test situation, it might be argued that the a priori model 
structure and parameter estimates may be an inaccurate description of the 
physical system to the extent that any claim to a precise control of the real 
output amplitudes is not valid. Still, a precise control on the output 
amplitudes for the a priori model provides a tool for conservative design of 
the input in terms of maximum output variable excursions, and thus gives 
the experiment designer the most control possible over the output variable 
amplitudes, given the a priori rhodel. With precise control over output 
amplitude excursions, any lack of confidence in an a priori model can be 
compensated for to a controllable and quantifiable extent by using 
conservative output amplitude constraints during the optimal input design. 

Since the input design for aircraft parameter estimation experiments 
must be done using a priori models and parameter values, it is in fact very 
likely that the a priori information does not precisely describe the physical 
situation (if this were not the case, motivation for the experiment would be 
lost). As a result, a flight test input design which is conservative as far as the 
output amplitudes produced, yet still effective relative to information content 
in the flight test data, is desirable. The task of keeping the aircraft response 
variables within limits is at odds with the objective of exciting aircraft motion 
as much as possible in order to obtain accurate parameter estimates. 

Including the output amplitude constraints directly in the input design 
problem formulation makes it possible to satisfy the output constraints while 
maximizing the information content in the data. This capability for 
conservative input design can help avoid ad hoc flight test procedure changes 
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during flight test operations, and also cut down on flight time requirements 
by assuring that data from flight test maneuvers can be analyzed using the 
same model structure used to design the input. For the case of multiple input 
design with higher order models, the direct inclusion of output constraints 
into the input design problem formulation is critical for the proper solution 
of the problem. The importance of output amplitude constraints to the 
proper design of the optimal input for parameter estimation experiments 
requires that these constraints be included as part of the problem formulation 
(equation (2-6)). 

In the present work, the admissible control set was comprised of square 
wave inputs. Square wave inputs are defined as controls whose amplitudes 
at any time are either full positive, full negative, or zero, so that the 
admissible control set, Uj, for the j lh component of the m-dimensional u 
vector is given by 

u j = { j / 0 / -Uj } (2-25) 


where |ij is the input amplitude constraint, and the j lh component of the 
u(t) vector, Uj(t), must take one of the values of Uj at any time, 

Uj(t) e Uj , j=l,2,...,m , te [0,T] q . 


With each of the m components of the control vector able to assume 
any of three distinct values, it follows that the admissible control set, U, for 
the m-dimensional control vector, u, contains 3 m distinct members. 


U =( wi,w 2 ,...,w 3 ") 


(2-27) 


where wj, W 2 , etc., are distinct m-dimensional control vectors. The control 
vector, u, must be equal to one of the members of the admissible control set 
for square wave inputs at any time, 

u(t) e U , t e [0,T] 
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The control vector u(t) represents perturbations about the input 
required for the nominal or trim condition at the start of the flight test. Thus 
the (generally nonzero) input associated with the nominal or trim condition 
is treated as zero deflection for the parameter estimation input. Limiting the 
possibilities for the input form is useful in the selection of the optimal input 
because it is easier to select the best input from the lot when the lot becomes 
smaller. The importance of this restriction on the input form will be made 
clearer in chapter III when the problem solution is described in detail. 

There are several reasons for the rather severe restriction described 
above concerning the allowable form of the optimal input for parameter 
estimation experiments. First, in cases where a human pilot implements the 
designed inputs, the square wave input was considered to have the highest 
repeatability and lowest potential for implementation distortion. This point 
is of course rendered moot in the event that the designed inputs are 
implemented by an automatic control system. Still, the simplicity of the 
square wave input form makes it easier to account for control system 
dynamics such as lag and hysteresis, using the solution method described in 
the next chapter. Such control system characteristics are present regardless of 
whether a machine or a human is at the controls. 

Careful reading of the references suggests the use of a simple square 
wave input for aircraft parameter estimation experiments. Chen [9] and 
Reid [3] showed that the minimum flight lest time for a simplified problem 
using a matrix norm of the information matrix as the cost function results in 
a "bang-bang" optimal control. Mulder [2] and Plaetschke and Schulz [7] 
evaluated several different input designs in flight, and found that a type of 
square wave input, the so-called "3211" input (see figure 2), was arguably the 
best performer overall, although this was not claimed by any of the authors. 
Gupta and Hall [4] used techniques developed by Mehra [6] to design an 
optimal input based on a scalar norm of the information matrix. As part of 
that study, the designed input was approximated by a square wave in order to 
make the input easier to implement. The authors compared the performance 
of the optimal input to its square wave approximation through simulation 
and found that the square wave approximation actually produced lower 
standard errors for some of the model parameters. Consideration of these 
results weighed heavily in the decision to restrict the input forms to square 
waves. 
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High information content in the data from a flight test is achieved by 
the use of abrupt, sharp-cornered inputs, since the frequency content for this 
type of input has a broad range, and thus is good for exciting the various 
modes of the aircraft response. These response modes have differing 
frequencies, which of course are not precisely known in the case of flight 
testing for the purpose of parameter estimation. Inputs comprised of 
frequencies covering a broad range are therefore superior to inputs with a 
more limited frequency content, such as a sinusoid. The input design, which 
necessarily is based on a priori models and parameter values, may not include 
exactly the proper frequencies to properly excite the response modes of the 
real physical system. Use of wide frequency band signals, such as square 
waves, is a small hedge against this problem. Interesting experimental data 
related to this issue has been presented by Mulder [2], 

Finally, the restriction on the form of the input is advantageous in the 
solution algorithm used to compute the optimal input in this study. This 
point will become clear in the next section, which is concerned with the 
solution method. 
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C hapter II I - S o lution Methodolog y 


The optimal input design problem, as formulated in the preceding 
chapter, was solved using the principles of dynamic programming [10] [11]. A 
brief general description of dynamic programming, as applied in this work, is 
included as Appendix B. In the form of dynamic programming used here, the 
outputs, sensitivities, and optimality criterion associated with each possible 
control at each time step are computed, followed by an efficient and 
sequential solution of the resulting high order combinatorial problem. In 
this way, the optimal input sequence, in terms of the lowest value of a chosen 
optimality criterion, can be found for any time period. 

There are considerations particular to optimal input design for aircraft 
parameter estimation which make a solution algorithm based on dynamic 
programming principles advantageous. These considerations will be 

highlighted in the course of the description of the solution algorithm, which 
follows. 

Equation (2-6) of Chapter II states that selected model output 
amplitudes are to be constrained within specified limits. For purposes of 
illustration, assume that only two such outputs, y, and y 2 , need to be 
constrained. The allowable output space at any given time then may be 
represented by a plane region whose borders correspond to the upper and 
lower amplitude constraints for y, and y 2 , see figure 3. The plane region is 
divided into discrete output space boxes. Time is also divided into discrete 
steps called stages. The time length of one stage is represented by T stage , and is 
assumed constant. The constrained outputs of the system are examined at 
every discrete time, separated by the time T stagc . Feasible continuous outputs 
at any time must be contained in one of the discrete output space boxes for 
that time. For the aircraft parameter estimation experiment, the initial 
condition of the constrained outputs is zero, since it is assumed that the 
model was derived relative to some specified condition. The initial output 
falls within a specific box in the discretized output space at the initial time. 

The fact that the initial condition for the problem is known means that only 
one sweep through time will be necessary for the solution. Assume that a 
scalar input is to be designed, so that the only input possibilities are full 
positive amplitude, full negative amplitude, and zero, for a square wave 
input form. Starting at the initial condition box in discretized output space. 
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all possible controls are applied over the time length T sta g e , and the 
consequences of each control possibility are computed. The consequences of 
each candidate control include the states, outputs, sensitivities, and cost 
associated with each control, deferring for the moment the question of how 
the cost is computed. Fourth order Runge-Kutta numerical integration of the 
equations of motion (2-1) and the state sensitivity equations (2-12) over the 
time step Tstage/ along with calculation of the outputs from equation (2-2) and 
the output sensitivities from equations (2-11), is required to determine the 
consequences of each control possibility. Over the period of time T stage / the 
input is constant and equal to one of the possibilities given by equations (2-25) 
through (2-28) with m=l. This amounts to allowing scalar square wave input 
sequences only. The resulting output values fall within one of the discrete 
output space boxes associated with the next stage, T stage later. The collection 
of all boxes which have been reached at the next time step might be called the 
reachable output space at that time. Using the state and state sensitivities of 
the system corresponding to each box in reachable output space at the next 
time stage as the new initial conditions, the process of computing the states, 
outputs, sensitivities, and cost associated with all possible controls is repeated 
to find the reachable output space for the following stage. If any of the 
possible controls applied over the time T s t a g e takes an output outside the 
feasible output space for the next stage, that control is excluded from 
consideration as part of the optimal input sequence. For figure 3, this means 
any control which produces an output beyond the output space plane region 
boundaries at any time would be discarded. This implements the output 
amplitude constraints in a simple and straightforward way. At this point it is 
easy to see the importance of keeping the number of possible controls to a 
minimum, since the consequences of each possible control must be computed 
for each box in reachable output space at each stage. 

For multiple input designs, there are simply more control input 
possibilities in the admissible control set, U, whose consequences must be 
computed for each box in reachable output space at each discrete time stage. 
For m inputs, there would be 3 m different control possibilities. The base three 
results from there being three choices for the value of each input vector 
component (full positive amplitude, full negative amplitude, and zero), and 
the exponent m is the number of elements in the control vector. 
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At each discrete time stage, the reachable output space is computed as 
the result of all possible controls starting at the reachable output space boxes 
found for the preceding time stage. After the computations proceed for 
several stages, some of the boxes in the reachable output space for a particular 
stage can be reached from more than one reachable output space box for the 
preceding stage, perhaps using different controls. Essentially, this means that 
more than one input sequence can be used to reach the same box in output 
space at a particular time. The discretization of the constrained output space 
was done to produce this occurrence. There must be some criterion for 
deciding which of the input sequences should be preferred to arrive at that 
box in output space at that time. This criterion is the minimization of a 
chosen function, which will be referred to as the value function, in keeping 
with standard dynamic programming terminology. Minimization of the 
value function is the optimality criterion for the optimal input design. 

For the minimum time problem, the objective of the input design is to 
reach specific goals for the Cramer-Rao bounds in the minimum time. 

Denote the goal for the Cramer-Rao bound of the k th parameter by , and 
the actual value of the Cramer-Rao bound obtained using a candidate input 
sequence up to the i‘ h time stage, by a kj for the k th parameter, k=l, 2, ..., p. 

Both o k . and £ k are positive values. Now define the value function as the 
square of the shortest Euclidean distance between the point in the p- 
dimensional space whose coordinates are the Cramer-Rao bounds at a given 
time, o k . , and the p-dimensional parallelepiped on and within which all the 

goals for Cramer-Rao bounds are achieved. The value function at the i th 
time stage, 'Fj(O) , for a candidate input sequence, fl, is given by 



where the depend on the a priori model and the candidate input sequence 

up to the i th time stage, and £ k are constant Cramer-Rao bound goals chosen 
a priori. 

Figure 4 depicts the value function calculation for the simple case of 
two model parameters (p=2). Cases A, B, and C in the figure represent three 
cases where a nonzero value function is possible for this simple situation. 
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The location of points A, B, and C would be determined by the computed 
Cramer-Rao bounds associated with, say, three different inputs, or perhaps 
the same input sequence at three different times. 

In order to use the value function given in (3-1), there must be some 
way to compute the Cramer-Rao bounds, a k . , sequentially. This is in contrast 

to the usual batch calculation of the dispersion matrix implied by equation 
(2-15). A sequential dispersion matrix, D|, calculated at the i th time stage, 
would give values as the square root of its diagonal elements (see equation 

(2-17)). D| would include information from a simulated experiment 
conducted with a candidate input sequence and the a priori model, up to the 
i lh time stage. 

The method for this sequential calculation is due to Chen [9], who 
applied the matrix inversion lemma [13] to equation (2-15) to produce a 
sequential update expression for the dispersion matrix in the form: 


^1+ 1 - D, - D, 


ay(i-^-l) \ T j 3y(i+l) \ D 

ae ) \ ae ) 1 


ay(i+l) 

ae 



ay(i+l) 

ae 



(3-2) 


Di+i is a symmetric matrix, since D must be symmetric (see equation (2-15)). 
Elements of the sensitivity matrix in (3-2) are obtained by solving the state 
sensitivity equations, (2-12), and the dynamic system equations, (2-1), and 
then using the output sensitivity equations, (2-11). 

The dispersion matrix update in (3-2) produces values of Okj +] f° r 

k=l,2,...,p, which give the coordinates of a point in p-dimensional space 
corresponding to a particular control sequence at the (i+l) th time stage. This 
sequential calculation allows the association of a scalar value function (from 
(3-1)) with each candidate input sequence in a manner which is stepwise in 
time. 

From numerical experimentation, it was found that the initial 
dispersion matrix, Di , could be made equal to a p-dimensional diagonal 
matrix having large diagonal elements, with no effect on the calculation of 
Dj+i past the first sequential update. Thus, the algorithm initializes Dj as 

Di = 10,000.0 (I p ) (3-3) 
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where I p is a pxp identity matrix. 

Each box in the reachable output space at each time stage has a value 
function, x Fj($) , associated with it, calculated from equation (3-1). As 
mentioned above, there may be more than one input sequence which results 
in output response values being contained within a single box in discretized 
output space. The input sequence which results in the lowest value of 'f'j(fl) 
is saved and associated with that particular box, along with the value of x Fj(d). 
Inferior input sequences with higher H'j(i3) values, which reach the same box 
in output space, are discarded. This sequential selective process is the 
manifestation of Bellman's principle of optimality [10] [11], and yields the 
global optimal input for any reachable output space box at any time stage in 
terms of the lowest value function. If each reachable discrete output space box 
is denoted by j, for j=l,2,...,n y (i), where n y (i) is the total number of reachable 
discrete output space boxes at time stage i, then the lowest value function to 
be associated with each reachable discrete output space box, j, at the i th time 
stage is denoted by and given by 


l Fi(i3*j) = min £ v Fi(Oj)J , j-l,2,...,ny(i) 

flje U 


(3-4) 


where the iDj input sequence up to the i*h stage produces output in the j 1 * 1 
reachable output space box, and the notation fyeU implies that the input 
sequence dj is comprised of elements of U which have remained constant for 
a period of at least one stage time of length T stag0 . The minimization in 
equation (3-4) is done by comparing all v F i (i') J ) for a given i and j, and selecting 
the minimum value. This minimization approach is feasible because the 
number of 'Fj(dj) for any given i and j is small, due to the discretization of the 
continuous physical problem. 

The fact that this solution process gives the input sequence with the 
global minimum value function for each reachable output space box at each 
time can be understood by realizing that this optimality is imposed 
sequentially from the initial time, and covers all possible (reachable) output 
space, subject to the chosen output space discretization. This latter condition 
means that all continuous outputs inside a discrete output space box are 
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considered the same point in output space. In this way, all of reachable 
output space at each time is associated with an input sequence $*j which 
reaches that output space at that time with a minimum value function, with 
the understanding that all continuous output space inside a discrete output 
space box is treated as the same output. 

For the minimum time problem, a zero value of corresponds to 

the achievement of all Cramer-Rao bound goals, since the k th term in the 
value function summation in equation (3-1) is included only when a k . > £k- 
This avoids imposing a penalty for o^. < Thus, to find the minimum time 

to attain specified Cramer-Rao bound goals, the calculations are stopped at the 
first time stage where any box, j, in reachable output space is associated with a 
zero value function, 'Pj(t'>*j). The first time stage with zero % FjO&*j) gives the 
value of the minimum time. The input sequence required to reach that box 
in reachable output space is the optimal control. Thus, for the minimum 
time problem, the cost, J, is 


J = T = (T stagl . ) i when a kj < Ck V k=l,2,...,p (3-5) 


where T s t age is the constant stage time. The minimum value of i which 
satisfies equation (3-5) is called n max , and this integer corresponds to the 
minimum test time required to achieve all Cramer-Rao bound goals. It 
follows from equations (3-1) and (3-5) that the integer n max is the index of the 
first stage which has H'jU'^j) = 0 for some je {l,2,...,n y (i)}. Since x Fi(0*j) should 
be a positive monotonically decreasing function as time goes on (and more 
information is added), the search for a zero value of 'Fi(f>*j) can be made by 
finding the smallest value of v Fi(ft*j) searching over all values of j at a given 
time stage i. This search will give the optimal value function attainable at 
any location in feasible output space at the i lh time stage. The optimal value 
function at stage i, %*, is given by 


min [^(OV] 

j e {l,2,...,n y (i)} 


(3-6) 
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So that 


n max = i when *¥* = 0 (3-7) 

The optical input sequence up to the i th stage is given by "0*j*(i) , where j*(i) 
is the index j of the optimal value function found from equation (3-6). Then 

Vi* = x Pi(d*j* {i )) (3-8) 

and for i=n max , 

X ^ n max = ^ n max^*j*( n max)) ( 3 ’ 9 ) 

with n max from equation (3-7) in conjunction with equation (3-6). 

For the fixed time problem, the solution method is the same, except 
that the target values for all the Cramer-Rao bounds are set to zero, i.e., 

Ck = 0 , k = l,2,...,p (3-10) 

and the cost is simply equal to the value function from equation (3-1), with all 
Ck - ke., 


J £ Kil 2 

k - 1 


(3-11) 


The optimal value function, H'j* , for each time stage, i, is determined in the 
same way as for the minimum time problem, that is, by equations (3-4) and 
(3-6). For the fixed time problem, the optimization is stopped after a fixed 
number of time stages, given by n ma x , where 

T 

n max “ -p 

st ag c (3-1 2 ) 


Here, T is the fixed test time chosen a priori, and T stagc is the constant value 
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of the stage time. The optimization is finished when the number of stages, i, 
reaches n ma x / computed from equation (3-12), 

T 

1 = n max = y 

sfagc (3-1 3) 

Equations (3-8) and (3-9) for the optimal input sequence hold exactly as before, 
with n max now given by equation (3-12). 

Equations (3-10) and (3-11) implement an implicit weighting on the 
parameter accuracies, in that the model parameters with larger absolute 
values tend to have larger Cramer-Rao bounds o k . associated with them, 

which means that their contribution to the cost in equation (3-11) is larger. In 
fact, any fixed test time formulation must address the problem of how to 
weight the relative importance of each parameter accuracy. The scheme used 
here was found to be most effective in the sense of achieving the lowest 
Cramer-Rao bounds, because although the optimizer pays more attention to 
reducing the larger Cramer-Rao bounds, these particular Cramer-Rao bounds 
also have the most potential for reduction, because the parameters associated 
with them are large. Large parameters in the model significantly affect the 
model output, which is equivalent to saying the sensitivity matrix elements 
associated with these parameters are large. The result is that the complete 
absence of weighting factors implied by equations (3-10) and (3-11) was found 
to be the most effective parameter weighting scheme for the fixed time 
problem. 

A fact glossed over until now is that in order to compute the chosen 
value function in (3-1) or (3-11) recursively, the state sensitivity equations 
(2-12) must be integrated numerically in time. These equations and the 
output sensitivity equations (2-11) require that the dynamic system equations 
(2-1) also be solved, since the state variables appear in the sensitivity 
equations. At each discrete time stage, a different input sequence, f)j*, will be 
associated with the j th reachable output space box. Thus each reachable box in 
output space has a different state and state sensitivity time history associated 
with it. In addition, the p(p+l)/2 unique elements of a sequentially updated 
symmetric dispersion matrix, Dj , will also be associated with each reachable 
box in output space, since the dispersion matrix is a function of the input 
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sequence through the sensitivities and the system dynamic equations. These 

considerations point toward a large memory requirement for the solution of 
the problem. 

However, there are aspects of the problem formulation which rescue 
what would otherwise be a bleak situation in terms of memory requirements. 
Since the dynamic system model is assumed to be derived relative to some 
specified condition, the initial conditions on the state and output variables 
are assumed known and constant for the experiment. The problem of 
designing inputs for a parameter estimation experiment differs from a 
general optimal control problem in that a "one time" input sequence (as 
opposed to a feedback control) is desired for only one initial condition. Thus, 
only one optimal input sequence beginning at the initial time must be saved 
for each reachable box in output space at each stage, and the propagation of 
the dynamic system equations, state sensitivity equations, and the unique 
elements of the dispersion matrix for each reachable box in output space can 
be done from stage to stage. Stated another way, the optimal input sequence 
for all time since the initial time must be saved for any particular reachable 
box in output space at any stage, but the time histories of the states, state 
sensitivities, and unique dispersion matrix elements can be discarded for all 
stages except the immediately preceding stage, since only these immediately 
preceding values are necessary to propagate the solution. The outputs and 
output sensitivities can be computed from the states and state sensitivities, 
and known quantities (see equations (2-2) and (2-11), respectively). Therefore, 
the variables that must be saved for each reachable box in output space for the 
current (i th ) stage only are given in the top part of Table 1. Since only a few 
input choices are allowed, each input choice can be associated with an integer 
index, so that the optimal input sequence, iy, saved for each reachable box in 
output space at any time stage is a string of integers, one integer for each stage 
from the initial time to the current time. This is shown at the bottom of 
Table 1, where array IOPTU holds the integer strings representing fy for each 
reachable output space box at each time stage. After the input design is 
completed, and the optimal input sequence determined, the optimal input 
sequence can be replayed using a numerical integration of the dynamic 
system equations along with the control system dynamics and the output 
equations, in order to reconstruct the state and output time histories 
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associated with the optimal input. This scheme keeps the memory 
requirement small. 

Several special features were incorporated into the basic algorithm 
described above in order to improve the computed solution. These features 
are outlined below. 

Selection of the output space discretization The chosen size of the output 
space boxes is important to the optimal input solution. The trade-off is one of 
increased memory and computation time requirements for more output 
space boxes versus a compromise of the optimal solution for too few boxes. 
The number of boxes to be used may be determined by solving the problem 
repeatedly with an increasingly finer grid for the constrained output 
variables. When the optimal input design is unchanged in going to a finer 
grid, the process may be stopped, and the last output space discretization used. 
The algorithm described in the present work has been designed to facilitate 
this type of calculation by allowing the number of divisions covering the 
amplitude range for each constrained output to be specified as input to the 
program. Combining this information with the output amplitude constraint 
values (rik in equation (2-6)), which must also be input to the program, 
completely specifies the discretization of the output space. If the number of 
discrete output space boxes for constrained output y k (t) is denoted by n yk , the 

length of the side of the discrete output space box associated with constrained 
output yk(t) is given by E k , 


n *k (3-14) 

Then, the dynamic programming grid for constrained outputs yk(t), 
ke {1,2,.. .,q), is set up by computing the boundary values for the discrete 
output space grid. These boundary values are stored in array YLIM, which is 
computed by 


YUM(k,j) = y k (0)-q k + (j-l)e k , j=l,2,...,(n yk +l) , ke {1,2, ...,q} (3-15) 
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where yk(0) is the known initial value of the k‘ h constrained output, and only 
the indices corresponding to constrained outputs are used for k. Equations 
(3-14) and (3-15) specify the dynamic programming grid boundaries for 
discrete constrained output space. 

S e l ection o f jh g - stage tim e and input frequency spectrum Since control can 
change only at the discrete time points separated by T slagC/ and only full 
positive, full negative, or zero input amplitudes are allowed on each control 
vector component, it follows that the available input frequency spectrum is 
influenced by the choice of T s t agc • In particular, if the high frequency 
components associated with the sharp corners of a square wave are 
disregarded, the larger the value of T sta g C , the more limited the input 
frequency spectrum palette. Very small values for T stagc give the freedom to 
change the control at more closely spaced times, due to the tight spacing of the 
stages. In this case, a wider, more nearly continuous frequency spectrum is 
available to the input design algorithm. A large T sta g C constrains control 
changes to be more widely separated in time, corresponding to lower, more 
widely spaced frequencies in the input. The trade-off is that smaller values of 
Tstagc mean more stages for a given time, which means a larger memory 
requirement, versus the better input frequency resolution available with 
smaller values of T stage . 

The algorithm automatically computes the stage time. Several issues 
must be weighed to choose the stage time appropriately. On one hand, a 
small stage time is desirable in order to admit a larger, more nearly 
continuous field of candidate input frequencies. However, limiting the stage 
time value on the low end is the memory capacity, particularly for large test 
times. Even more limiting on the low end for the value of the stage time is 
what might be called discretization error. If the stage time is very small for a 
fixed control amplitude, it may happen that application of all control 
possibilities over one stage time produces output which remains inside the 
originating box in output space, regardless of the control applied. This 
happens most apparently at the outset of a single input design, when full 
positive and full negative controls produce equal value functions. In effect, 
the algorithm would discard the second of the applied controls (full positive 
or full negative) as inferior, and in so doing would discard half the 
subsequent decision tree, and thus seriously compromise the optimal 


33 



solution. It is therefore essential that at least one nonzero control application 
at the outset produce output which travels outside the originating box in 
output space, and this criterion is used by the algorithm to determine the 
stage time. For the chosen output space discretization, the equations of 
motion are integrated from the initial conditions until at least one output 
crosses a box boundary in discretized output space. The time for this 
occurrence is designated as the stage time, and calculations proceed as usual. 
Assume that the initial value of the k th constrained output is between 
YLIM(k,jL) and YLIM(k,jH), that is, 

YLIM(k,jL)<y k (0)<YLIM(k,j H ) , ke {l,2,...,q} (3-16) 

Now, the dynamic system equations (2-1) are integrated forward in time using 
the first nonzero input from equation (2-27), and the resulting outputs are 
computed from (2-2) until 

y k (t) < YLIM(k,jL) OR y k (t) > YLIM(k,j H ) , ke {1,2,.. „q) (3-17) 

The minimum time for the condition in equation (3-17) to be satisfied for any 
constrained output, k, is designated the stage time, T s t a g C . So that 

T s tagc - minimum t satisfying (3-17) (3-18) 

With this feature, the minimum stage time possible (equivalently, the widest 
and most nearly continuous input frequency spectrum) is used for a given 
selected output space discretization. 

Practical upper limit on input frequency If the optimal input design is to be 
implemented by a pilot, an upper limit exists for the input frequencies which 
can be implemented, due to human physiology. The algorithm allows 
specification of this upper limit on frequency in the time domain by 
providing the means to specify the minimum time for the spacing between 
control changes. This time is an input to the program, and is denoted t mp . 
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The program converts this time into an equivalent number of stages, n mp , by 
the relation 


“mp t T 4 

stage (3- 19 ) 

where n m p is an integer, and one count is added to offset the truncation which 
occurs when assigning a real quotient to an integer variable. The top of Table 
1 shows that array IREPU keeps track of important control parameters, 
including the number of times the last control was applied consecutively for 
each reachable output space box at the current time stage. To implement the 
high frequency limit in the time domain, the program simply limits the 
admissible control set for n m p stages after a control switch to only the last 
control, i.e.. 


U L( j) = { w L{j) ) , L(j)e {1,2,.. .,3 m } (3-20) 

where L(j) is the integer index for the last control used to reach the discrete 
output space box, j, at the current time stage, which is the initial condition for 
the trial of control possibilities. If n^) represents the number of stages control 
w L (j) has been applied consecutively up to the current stage for the j th 
discrete output space box, the admissible control set is given by (3-20) as long 
as 


n L(j) < r» mp , je{l,2,...,n y (i)} (3-21) 

where i is the stage index, and n y (i) is the number of reachable output space 
boxes for the i^ stage. When equation (3-21) no longer holds, the admissible 
control set reverts to the usual set given in equation (2-27). In this way, an 
optimal input can be designed including the practical constraint imposed by 
the speed of reliable human pilot inputs. A similar limitation may be desired 
when using automatic implementation of the optimal input design for a 
variety of other reasons, such as instrumentation dynamic response, known 
control system dynamics, modelling errors related to frequency, or avoiding 
resonance with aircraft structural modes. 
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Control system dynamics A basic objection to a true square wave input 
design is that although the simplicity of its implementation is often touted, in 
practice, the implementation of a true square wave is impossible, due to the 
requirement for an infinite control surface displacement rate, see figures 1 
and 2. In order to design a truly practical optimal input, the capacity to 
include control system dynamics was built into the algorithm. The examples 
in chapter V employ a simple first order lag on the control surface deflections, 
though other more complicated dynamics could be incorporated. Referring 
again to the top of Table 1, the input amplitudes resulting from the 
application of the input sequence f)j* , including control system dynamics, are 
saved for each reachable output space box of the current stage. These 
amplitudes represent the initial control amplitudes for the computations 
over the next time period of length T stage . The control amplitudes 
represented by the components of each member of U in equation (2-27) can be 
thought of as commanded values for the control amplitudes. A scalar 
function of time, called AMPI.AG(t), represents the control system dynamics. 
Assuming wj is the control possibility being examined, and denoting the 
saved input amplitude vector for the current stage from array UASTR (see 
Table 1) by u 0 id, the input vector, u(t), is given by 


u(t) = (wi - u„id )AMPLAG(t - (i-l)Tstagc) + u„id (3-22) 


For a first order lag with time constant x, the AMPLAG(t) function is 

AMPLAG(t) = 1 - e (3-23) 

The control surface dynamics are included directly in the calculation of the 
value function associated with each control applied over any stage, and are 
not applied after the fact to make the implementation of the input design 
practical. This procedure preserves the optimality of the designed input; the 
input simply has added practical constraints imposed by the dynamics of the 
control system. 

Specification of initial and final input amplitudes As mentioned previously, 
the input design for aircraft parameter estimation experiments is essentially 
an input time history superimposed on trim or nominal input amplitudes. 
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Thus, a desirable feature of an input time history designed for aircraft 
parameter estimation is that the initial and final input amplitude values be 
zero. This enhances the safety of the flight test by helping the pilot recover 
the aircraft from the maneuver. It makes sense to include the requirement 
that the parameter estimation input begin and end with zero amplitude as an 
integral part of the optimal input design problem formulation. The 
algorithm in the present work is capable of implementing this requirement. 
The forward time evolution is continued until there is a box in output space 
which satisfies the conditions of optimality given previously, but also has 
been reached by application of the zero control vector. The time that the final 
control should be zero is an input to the program, called tf c . This time is 
converted to a number of stages, nf c , by 

n fc = -~ — + 1 

T smgc ( 3 - 24 ) 

where the same comments apply here as to equation (3-19) above. Now in 
addition to the stopping criterion for the minimum time problem (equation 
(3-7)) and the fixed time problem (equation (3-13)), another criterion is added 
for an optimal input ending with zero control deflections. Denote the index 
for the zero control vector in U by lo , l(>e {1,2,.. .,3 m }, and the control index of 
the last control for the j th discrete output space box at the current stage by L(j). 
The number of times this last control was applied is nL(j>. The quantities L(j) 
and n^) are saved for the current stage (see top of Table 1). Then the 
additional stopping criterion is 

L(j) = lo AND n L (j) > n (c , je {l,2,...,n y (i)) (3-25) 

where i is the stage index, and n y (i) is the number of reachable output space 
boxes for the i th stage. The first stage where this modified stopping criterion 
is satisfied is the end of the optimization, and all else remains the same. Due 
to the included control system dynamics, several time stages with a zero 
control may be required for the input amplitudes to actually return to zero. 

All the examples presented in the sequel include this requirement as part of 
the optimal input design. 
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Automatic adjustment of input amplitude Appropriate values for output 
amplitude constraints are easily specified a priori from physical 
considerations related to the safety of the aircraft and pilot, and the validity of 
the model. Input amplitude constraint values are not as easily chosen, since 
they may be strongly related to the imposed output amplitude constraints. 

For the aircraft parameter estimation problem, the starting point for input 
amplitude constraint values are the physical limits on control deflection, as 
implemented by mechanical stops, flight control software limiters, or linear 
control effectiveness. In general, higher input amplitudes give more system 
excitation per unit time, so that the accuracy requirements for the parameter 
estimates are achieved in a shorter time or more parameter accuracy is 
available for a fixed time. Thus, if no output amplitude constraints were 
imposed, the maximum input amplitude possible would be preferred. When 
output constraints are imposed, as is commonly the case, it may occur that the 
input amplitude is so large that feasible outputs are produced only when the 
controls are applied at the initial time. For example, assume that all controls 
applied at the second stage, starting from an output space box different from 
the initial box, produce output which travels outside feasible output space. In 
this case, all controls are disallowed, there are no reachable output space boxes 
for the next stage, and no optimal control exists. This problem is exacerbated 
by the maximum input frequency limit feature described above, since that 
feature specifies a minimum pulse width for each control. The algorithm has 
the capability to recognize this problem. This is done by setting a flag when 
the only feasible outputs computed for any stage beyond the second stage 
originate from the initial output condition. Corrective action consists of 
decreasing the input amplitude constraint by an amount which can be 
specified, and restarting the algorithm. Denote the input amplitude 
decrements by 8|ij for j=l,2,...,m. These input amplitude decrements are 
positive-valued inputs to the program. Then the input amplitude constraints 
are modified by 


Mj = Mj - 8 Mj 


j=l,2,...,m (3-26) 
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At this point, the algorithm is restarted. This iterative procedure ensures that 
the maximum possible input amplitude is used for the input design, 
consistent with the imposed output amplitude constraints. The input 
amplitudes are determined directly and iteratively as a function of imposed 
constraints on the output amplitudes and the practical limits for the 
maximum input amplitudes, all of these being easily specified in practice. 
Another occurrence which requires equation (3-26) can be called a monotonic 
optimal value function violation. This situation arises when the optimal 
value function for the i^ stage, H / j from equation (3-6), is substantially larger 
than 'Fj-i , the optimal value function from the preceding stage. In this case, 
the imposed output amplitude constraints are severely compromising the 
optimal solution. The problem is also addressed by reducing the input 
amplitude constraint using equation (3-26), and restarting the algorithm. 

Multiple input switching Multiple input design for aircraft parameter 
estimation experiments presents another practical difficulty, namely the 
coordination of the controls by a human pilot. In particular, a pilot may be 
unable to implement an input design where more than one input is changing 
at the same time in a relatively complicated fashion. To combat this problem, 
the optimal input design algorithm provides the capability to specify time 
periods where only one control will be allowed to vary from zero. This 
amounts to additional restrictions on the allowable controls at any given time 
stage and is an additional constraint on the form of the input. Assume that 
the control vector, u, has only two components (m=2). The time length 
during which only a single control vector component is allowed to move is 
an input to the program, called t sw . This time is converted to an equivalent 
number of stages in the manner described previously. 


n 


SW 


^sw 

T 

stage 


+ 1 


(3-27) 


Now denote the subset of the admissible control set U which moves only one 
control vector component, including the zero control vector, by Uj. Then, 


UdU| 


(3-28) 
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where in the case under consideration (m=2), U has 3 2 elements, and 



-Hi 

0 


(3-29) 


Ui represents the admissible control set for square wave inputs moving only 
the first component of the control vector. As long as 


i £ n sw (3-30) 

where i is the stage counter, then the admissible control set is given by Uj. 
When i > n sw , the admissible control set changes to U 2 , where 

U 2 

The extension to more than two components in the control vector is 
straightforward. The control switching time length, t sw , is constant for all 
control vector components, and the control sequencing is fixed, in the interest 
of simplicity. Modifications could easily be made in the direction of more 
sophistication. This additional constraint on the form of the optimal input 
(or any additional input form constraint, for that matter) actually decreases 
the computation time for the optimal input design, since fewer control 
choices must be considered for any reachable output space box at any time 
stage. 
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The following list with brief descriptions summarizes the features of 
the algorithm developed in this work for designing optimal inputs for aircraft 
parameter estimation. 

Global optimum solution - Expensive flight test time is minimized by 
formulating the problem as a minimum time problem. The global 
minimum time solution is the input sequence for the shortest flight test 
required to reach target values for the Cramer-Rao bounds. The global 
optimum fixed time input design is achieved by setting target values for the 
Cramer-Rao bounds to zero, and choosing the input sequence associated with 
the lowest Cramer-Rao bounds at the final time. Both solutions are globally 
optimal by virtue of Bellman's principle of optimality [10] [11], which is the 
modus opevundi of the algorithm. As with any other optimal solution, the 
optimality is meaningful only in the context of the problem formulation, 
which was given in Chapter II. 

Multiple input desi gn capability - Multiple input problems increase the 
number of input choices whose consequences must be evaluated for each 
reachable output space box at each time stage. The result is an increase in the 
computation time required to solve the optimal input problem. This is to be 
expected, since the multiple input problem is inherently more complicated. 
The solution algorithm has been formulated so that the additional memory 
required to solve a multiple input problem is small. 

Closed loop model capability - The program was designed in a modular 
fashion, so that various linear system dynamic and output models could be 
used. The algorithm can accommodate a model structure for linear closed 
loop models, thus allowing the design of optimal inputs for open loop 
airframe parameters when only closed loop data are available. 

Practically achievable input design - Constraints on the input form were 
incorporated to alleviate repeatability and coordination problems which 
might be encountered when a human pilot is used to implement flight test 
inputs, particularly for a multiple input design. In addition, the algorithm 
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incudes control system dynamics directly in the calculation of the optimal 
input design. This addresses the problem of implementation distortion, 
which would otherwise degrade the efficacy of the designed input. 

Direct implementation of output amplitude constraints - Inclusion of the 
output amplitude constraints in the problem formulation results in an input 
design which excites the system response modes to the fullest extent possible, 
but no further. This characteristic is important for a optimal input design in 
the practical sense. 

Input frequency chosen in the time domain - Control changes may occur only 
at discrete time points which are separated by one stage time. For square 
wave input forms, this means that the input frequency spectrum will be 
chosen by the determination of the optimal time for the input to switch. 
These switching times for the input will be optimal to within plus or minus 
one stage time. The algorithm automatically determines the minimum 
usable stage time, based on the output space discretization. This ensures that 
the widest and most nearly continuous input frequency spectrum is 
considered in the design of the optimal input.. Practical co ns traints on the a 
input frequency spectrum resulting from inherent characteristics of the pilot, 
flight test operational procedure, and the control system dynamics are 
implemented in the time domain. 

Single pass solution - The solution algorithm described here requires no 
iterative calculations or special start-up procedures. The optimal input 
solution is obtained in a single pass, which marches stepwise forward in time. 
When the problem solution is not possible with the specified input and 
output constraints, the algorithm decreases the input amplitude and restarts 
the solution. However, this iteration is on the proper conditions for the 
problem solution, and not on the solution per se. 

Well-suited to design studies - The algorithm presented here can be used for a 
meaningful trade-off design study of optimal input design for aircraft 
parameter estimation. For example, the analyst can easily examine the effect 
of changes in measurement noise characteristics, Cramer-Rao bound goals, 
input amplitude constraints, output amplitude constraints, or control system 
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dynamics in terms of their effect on the required flight test time or achievable 
parameter accuracies. 

To complete the picture, drawbacks of the algorithm are listed and 
described below. 

M emory requirement - A requirement for a large amount of memory is 
inherent in dynamic programming solutions to optimal control problems. 
The originator of the dynamic programming approach, Richard Bellman, has 
referred to this problem as " the curse of dimensionality " [10] [11]. In the 
present work, the memory requirements were kept to reasonable levels by 
taking advantage of characteristics peculiar to the problem of optimal input 
design for aircraft parameter estimation. These characteristics are the "one 
time nature of the input design, the fact that the dynamic model was 
formulated relative to a specified initial condition, and the limitations which 
could be legitimately imposed on the form of the input. For this work, 
memory requirements were impacted most significantly by the number of 
constrained output variables required, the output space discretization, and the 
maximum number of time stages. To a much lesser extent, the model order 
and the number of parameters in the model increased memory requirements 
as well. While more complex problems will naturally increase memory 
requirements, this is not seen as a serious problem for several reasons. First, 
the algorithm has been designed with an eye toward minimizing the memory 
requirements, in a time when the run-time memory of new computers is 
increasing at a rapid pace. In addition, the current version of the program can 
handle a sixth order dynamic system model, with twelve model parameters, 
two inputs, two constrained outputs, and a maximum of one hundred fifty 
time stages, with still some room for expansion. The current capability will 
be sufficient for many practical problems of interest. 


Ru n tim<? requirement - The required run time depends mainly on the model 
order, the number of model parameters, the output space discretization, and 
the number of control choices available for each output space box at each time 
stage. To give some idea of the execution time required, using a Digital 
Vaxstation II, the open loop single input example presented in this work took 
about 3 minutes for a 3.2 second input design. 


43 




The first example was investigated previously by Mehra [5] and 
subsequently by Chen [9]. The problem is to design a single input for the 
purpose of estimating the parameters contained in an approximate model of 
the short period dynamics for an airplane. The model is given by 



where a is the angle of attack, q is the pitch rate, and 5 C is the elevator 
deflection. The output model is . 

I _ f 1 0 l|ct(t)| 

y 2 (t)J 1° lJlq(t)J 


The measurement equations are 
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where u(i) is the i lh realization of a zero-mean Gaussian white noise vector 
random process with the following measurement noise covariance matrix: 
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The a priori values of the model parameters are the true values for the 
parameters, given in Table 2. The physical system is the system given above 
using the true parameter values from Table 2, which means that there is no 
modelling error in the a priori model. For this problem, Mehra [5] used a 
fixed time for the flight test data run, and maximized the trace of the 
information matrix ( see equation (2-9) ), subject to an energy constraint on 
the input of the form given in equation (2-21), i.e.: 



(5-6) 


where T=4.0 seconds and E=311. Figure 5 shows the Mehra input design, and 
the associated time histories for a and q. The third column of Table 2 gives 
the performance of this input design in terms of the Cramer-Rao bounds. 

In reference [9], Chen solved the identical problem using the same 
energy constraint on the input. Instead of minimizing some scalar norm of 
the information matrix, the minimum time required to achieve the 
parameter accuracy represented by the Cramer-Rao bounds from the Mehra 
input was sought. The procedure used was iterative, and assumed the input 
form to be a member of an orthogonal set of Walsh functions, which are 
"bang-bang" input forms. The solution obtained was suboptimal, and is 
shown in figure 6, with resulting Cramer-Rao bounds given in the fourth 
column of Table 2. 

The new technique for optimal input design was used to find the 
minimum time solution for square wave input forms for the same problem. 
Since no energy constraint of the form (5-6) was imposed, the problem was 
solved for two cases - once using the maximum input amplitude of the Chen 
solution as an input amplitude constraint, and a second time using the 
maximum input amplitude of the Mehra solution. The results for the 
former case are shown in figure 7, while the results for the latter case are 
given in figure 8. Values for the Cramer-Rao bounds for these inputs appear 
in columns 5 and 6 of Table 2. 

The problem was solved a third time using the technique developed in 
the present work. Input amplitude was again constrained to be the same as 
the Mehra solution, but this time a first order control system lag was 
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included, with time constant equal to 0.1 second. The input design and 
output variable time histories are given in figure 9. The resulting Cramer- 
Rao bounds for the parameters are given in column 7 of Table 2. 

The Cramer-Rao bounds in Table 2 for the input designs of Mehra and 
Chen were obtained with the same computing means used to obtain the 
Cramer-Rao bounds for the optimal input designs generated by the technique 
described in this dissertation. That is, only the form of the input was 
borrowed from the works by Mehra and Chen ; the resulting Cramer-Rao 
bounds were computed with the same computer and algorithm ( equation te- 
ll) and numerical solution of (2-1) and (2-12) ). This was done in the interest 
of a fair comparison among the input designs. 

The input designs associated with the present work included the 
stipulations that the input begin and end with zero amplitude, and the 
minimum pulse width was set at 0.6 second. In addition, the ot and q output 
variables were constrained to ± 10 degrees and ±12 degrees per second 
respectively, so as to limit output variable excursions to values similar to 
those obtained from the Mehra and Chen input designs. The sampling 
interval. At, was set at 0.02 second. The goals for the Cramer-Rao bounds 
were taken to be the values associated with the Mehra input, from column 3 
of Table 2. Thus, the objective for the optimal input designs was to attain 
these goals in minimum time. 

Column 5 of Table 2 shows the results of the present work using the 
same input amplitude as for the Chen design. A lower value of total time is 
attained by giving up some accuracy on the Z force parameters, while picking 
up some additional accuracy on the moment equation parameters. All 
Cramer-Rao bounds are lower than the Mehra input values, as required. This 
trade-off among the parameter accuracies beneath the goal values in order to 
achieve a smaller test time was done implicitly and automatically by the 
algorithm, thus obviating the need for parameter weighting and all the ad 
hoc procedures that usually accompany it. Maximum output amplitudes 
were virtually the same as for the Chen input design, and lower than those 
for the Mehra input design, as required by the output constraints built into 
the solution algorithm (see figures 5,6, and 7). 

The results in column 6 of Table 2 demonstrate that if the higher 
maximum input amplitude of the Mehra design is used as an input 
amplitude constraint, the goals for the Cramer-Rao bounds can be achieved in 
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an even smaller total time. All Cramer-Rao bounds were below their 
respective goal values from the Mehra input in column 3 of Table 2. Output 
variable amplitudes were again within their constrained ranges, see figure 8. 
This case demonstrates some of the points made in chapter II regarding input 
characteristics and the appropriate types of constraints for aircraft parameter 
estimation experiments. 

The effect of including a control system lag on the input design 
solution is given in the last column of Table 2. Cramer-Rao bounds were 
again less than the goal values from the Mehra input for every parameter, as 
required. These accuracies were virtually the same as in column 6, but 
required a test time increase of 0.16 second. This slightly longer test time was 
expected, since the higher frequencies associated with the sharp corners of a 
square wave input were filtered out by the control system lag. Output 
variable amplitudes were again within their constrained ranges, see figure 9. 

Table 3 shows the results of maximum likelihood estimation of the 
model parameters, based on simulated flight test runs using each of the input 
designs discussed above and listed in Table 2. Excellent agreement was 
obtained between the parameter standard errors in Table 3, and the Cramer- 
Rao lower bounds for these standard errors, given in Table 2. This 
exceptional agreement was due mainly to the fact that no modelling error 
existed for the simulated flight tests and the measurement noise was precisely 
known. Real flight test situations do not enjoy such luxuries. As a result, the 
performance of the input designs would be degraded in a real flight test 
situation. However, Mulder [2] has presented experimental evidence that the 
relative merits of various input designs based on the Cramer-Rao bounds 
remain intact when applied in a real flight test environment. This was 
demonstrated likewise for the simulated flight test results in Table 3. 

The results in Tables 2 and 3 and figures 5 through 9 demonstrate the 
quality of the input designs which result from the problem formulation and 
solution algorithm described in this dissertation. 
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Example 2 


This example exhibits the expanded capability of the input design 
algorithm by computing the optimal input for a fourth-order model with 
twelve model parameters ancT two inputs, subject to restrictive output 
amplitude constraints. 

The lateral dynamics of an advanced fighter aircraft in level flight at 
10,000 m altitude and an airspeed of 179.7 mps may be represented by the 
following dynamic system and measurement model : 
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where P is the sideslip angle, p is the roll rate, r is the yaw rate, <f> is the roll 
angle, 8 a is the aileron deflection, and 8 r is the rudder deflection. The output 
is given by r i 
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The constraints are 




0.07 radians V t 
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< 0.15 radians V t 
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The measurement equations are 


’y mi (o' 

y m2 (i > 

y m3 (i) 

y m4 (i) 


10 0 0 
0 10 0 
0 0 10 
0 0 0 1 


>!«)■ 
yP + 

y 3 (t) 

y 4 (0 


'i)i(i)" 

u 2 (i) 

t)3(i) 

t) 4 (i) 

- 




(5-14) 


where \)(i) is the i 1 ^ realization of a zero-mean Gaussian white noise vector 
random process with the following measurement noise covariance matrix: 
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The values for the diagonal elements of R in equation (5-15) are derived from 
information in reference [14] concerning appropriate data system 
characteristics for aircraft parameter estimation experiments. These values 
represent random fluctuations due to the instrumentation system only, and 
thus do not include such real effects as modelling error and gusts, which 
typically manifest themselves in the measurement noise. 

As in the last example, the physical system is represented by the model 
structure given above using the true values of the parameters, which are also 
assumed to be the a priori parameter values. The true values of the model 
parameters are given in Table 4. The above model does not correspond to any 
existing fighter aircraft or instrumentation system. The nominal values used 
are for demonstration purposes only, although the model structure is quite 
generally applicable. 
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The output amplitude constraints imposed on the sideslip angle and 
the roll angle are such that a straightforward input design using past 
experience and heuristic arguments is difficult. A common practical input 
design procedure employs doublets in the rudder and aileron with 
frequencies close to the damped natural frequency of the aircraft Dutch roll 
mode. This type of input is shown in figure 10. A first order lag with time 
constant 0.1 second was implemented to model pilot and control system 
actuation delay. By trial and error, an input amplitude of 0.07 radian was 
found to satisfy the output amplitude constraints using the doublet inputs. 
Cramer-Rao bounds associated with this 10 second input design appear in 
column 3 of Table 4. 

For the optimal input designs, input amplitude constraints were 
applied instead of input energy constraints, as discussed previously. The 
same input amplitude constraints were used for both control surfaces (see 
equations (5-9) and (5-10)), and the controls were sequenced, with the control 
switching time, t sw , set to 5 seconds, and the rudder control sequenced first. 
Maximum input frequency depends on the minimum pulse width. For this 
problem, the assumption was that 0.6 seconds was the minimum time 
necessary to acquire a control amplitude and return the control to zero (i.e., 
nominal or trim) deflection. 

In order to allow comparison with the doublet inputs, the optimal 
input design technique was used to produce a fixed time design by setting the 
goals for all Cramer-Rao bounds to zero, and limiting the run of the input 
design program to a maximum of 10 seconds. The resulting input design and 
output responses appear in figure 1 1 . The Cramer-Rao bounds were lowered 
for all parameters (compare columns 3 and 4 of Table 4), with all output 
amplitude constraints satisfied. 

In figure 12, the optimal input design was performed for the same 
problem using a larger (0.1 radian) input amplitude constraint on both inputs. 
Allowing a larger input amplitude constraint means that there is potential for 
more system excitation over a fixed time, compared to a case with a lower 
input amplitude constraint. Thus, one would expect a shorter test time to be 
required to reach the same Cramer-Rao bound goals, provided that the output 
amplitude constraints can be enforced with the larger input amplitude. For 
the case shown here, a minimum time solution was computed. The goals for 
the Cramer-Rao bounds were the accuracies associated with the doublet 
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inputs in column 3 of Table 4. The same accuracy or better for each parameter 
relative to the doublet inputs case was obtained in a shorter total test time of 
8.7 seconds, as shown in column 5 of Table 4. 

For all input designs presented thus far, the control inputs were 
sequenced; that is, only one control was moved at a time. This feature is 
helpful when human pilots must implement the input design. For cases 
when a computer can realize the inputs, the requirement for sequenced 
control inputs can be relaxed. This case is shown in figure 13, using the same 
input amplitude as in figure 12. Again, a fixed time solution was computed 
by setting all goals for the Cramer-Rao bounds to zero with a 10 second 
maximum test time. Columns 5 and 6 of Table 4 show that a gain of at least 
25% in accuracy for all parameters except those associated with the Y force was 
achieved for a 10 second total test time by relaxing the control sequencing 
requirement. The improvement in the Cramer-Rao bounds for column 6 of 
Table 4 relative to the Cramer-Rao bounds for the doublet inputs ranges from 
19.0% for Yp to 72.4% for Np, based on the Cramer-Rao bounds for the doublet 
inputs. All optimal input designs excited the aircraft response as much as 
possible to obtain small Cramer-Rao bounds, but still kept the output 
response within the imposed amplitude constraints. 

This example demonstrates the computation of practical, optimal, 
multiple input designs. In addition, the effects of increased input amplitude 
constraint values and relaxing the input sequencing were easily quantified in 
terms of lower Cramer-Rao bounds for a fixed time design, or reduced flight 
test time required to achieve fixed target values for the Cramer-Rao bounds 
for a minimum time design. The effects were precisely quantifiable because 
of the global optimality of the solutions. 
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Exam ple 3 


In this final example, the optimal input design technique is used under 
conditions which mimic a real flight test situation. The test aircraft is a six 
degree of freedom nonlinear simulation of the F-18 fighter aircraft [15]. Only 
the longitudinal dynamics of the aircraft are considered. 

The F-18 aircraft is designed with inherent longitudinal static 
instability, in order to enhance maneuverability and performance^ As a result 
of this, the aircraft must operate under closed loop automatic control at all 
times during the flight test for safety considerations. 

The flight condition is 20,630 ft. altitude and 8 degrees angle of attack at 
trim. Geometry and mass characteristics are given in Table 5. The open loop 
dynamic model is: 
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where V is the airspeed, a is the angle of attack, q is the pitch rate, and 5h is 
the symmetric stabilator deflection. The output is 
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The measurement equations are 
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where v(i) is the realization of a zero-mean Gaussian white noise vector 
random process with the following measurement noise covariance matrix: 
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A problem in designing an optimal input for this open loop model is 
that 5h , the stabilator deflection, is not under direct control of the designer. 
The deflection of this control surface is a function of flight condition and 
motion variables, as well as the pilot stick input, due to the action of the 
control system. It is therefore necessary to model the control system, so that 
the input design algorithm can deflect the stabilator indirectly in a manner 
which will be most advantageous for accurate estimation of the open loop 
parameters in equation (5-16). An adequate model structure for the control 
system was found to be: 
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where ri e is the longitudinal stick deflection at the pilot station. The output is 


y 4 (t) = 5h(t) 


(5-24) 


The measurement equation is 


where 


y m4 (i) = y 4 (i) + u 4 (i) i = 1, 2, ..., N 
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For the control system model, longitudinal stick deflection, r| e , as well 
as the aircraft motion variables V, a, and q, are treated as known inputs. This 
simple first order model for the control system is adequate for the small 
perturbations in T| e from trim which are typical of inputs designed for 
parameter estimation experiments. 

In general, the F-18 control system schedules leading edge and trailing 
edge flap deflections with angle of attack. This means that there should be 
two more rows in the control system model of equation (5-23) and two more 
controls in the open loop model of equation (5-16), corresponding to leading 
edge and trailing edge flap deflections. Including these controls, however, 
produces a collinearity in the data, since the flap deflections depend linearly 
on angle of attack at this flight condition. This causes parameter 
identifiability problems which are difficult to surmount. If one does not 
account for the fact that the flap settings change with angle of attack, a 
parameter estimation algorithm such as maximum likelihood assigns the 
effect of the flaps to the angle of attack parameters, thus rendering these 
parameter values inaccurate. For these reasons, the flap scheduling was 
turned off by setting the computed flap deflections from the flight control 
system subroutine inside the F-18 simulation to zero. With this done, the 
control system deflected only the symmetric stabilator, 5h , in response to 
longitudinal stick deflections. 

Combining the open loop dynamic model in equation (5-16) with the 
control system model in equation (5-23), and assuming that the measurement 
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noise in (5-20) is uncorrelated with the measurement noise in (5-25), gives the 
model for the closed loop longitudinal dynamics. The stabilator deflection is 
treated as an element of an augmented state vector, and can be viewed as an 
additional state, since its dynamics depend on state variables and the input, 
which is longitudinal stick deflection, rj e . In the notation already introduced, 
the closed loop model takes the form: 
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subject to the constraints 
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The measurement equations are 
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where 
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The open loop dynamic model (5-16) and the control system model 
(5-23) must be considered separately to estimate their respective parameters 
using maximum likelihood estimation. This is because the maximum 
likelihood parameter estimation algorithm integrates the dynamic and 
sensitivity equations using candidate model parameters as part of the process 
of finding parameter estimates. When attempting to estimate all the 
parameters in the closed loop model (5-27) using the data from a simulated 
experiment, the maximum likelihood algorithm diverges for nearly any 
initial values of the parameters, unless these initial values are very close to 
the "true" parameter values. This happens because the open loop system is 
inherently unstable. Estimating parameters in the open loop model (5-16) 
treating the closed loop stabilator from the experiment as a known input, and 
then separately estimating the parameters in the control system model (5-23), 
treating T| c , V, a, and q from the experiment as known inputs, eliminates the 
problem of a divergent parameter estimation algorithm. This separated 
estimation procedure was used to obtain all parameter estimates from 
simulated experimental data in this example. 

In accord with what might be true for a real flight test input design 
problem, it was assumed that the only available data from which to assemble 
an a priori model for the input design were from a simple doublet input 
sequence for the longitudinal stick deflection, q c . This doublet was in fact 
taken from an actual pilot input during the flight test of the real F-18 aircraft, 
see Klein et al [16]. The simple doublet input and the output time histories 
which result from applying this input form to the F-18 simulation test aircraft 
appear as figure 14. A first order lag with time constant 0.05 second was used 
in producing this input form. The same control input dynamics were used 
for all input forms in this example. 
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Based on the information from the simple doublet input, parameter 
estimates for the closed loop model were determined using maximum 
likelihood estimation and the separated model procedure described above. 
The dotted lines in figures 14(b) through 14(e) indicate model fits for each 
variable, as computed using the separated model procedure. The parameter 
estimates from the doublet input are given in Table 6. These parameter 
estimates and the model structure in equation (5-27) comprised the a priori 
model for the optimal input design. Figures 14(0 through 14(i) show the 
model responses using the closed loop model (5-27) and the parameter values 
from Table 6, as compared to the open loop model responses, which are the 
same as the dotted lines in figures 14(b) through 14(e). 

The parameters of principal interest in the open loop dynamic model 
(5-16) are contained in the bottom two rows, and are generally associated 
with the two state approximation of the short period dynamics. For the 
doublet input, the drag equation (top row) of (5-16) is included in the open 
loop model for two reasons. First, it was discovered that the control system 
model fit and predictive capability were considerably improved by including 
the velocity term in the control system model (5-23). This control system 
model fit is shown in figure 14(b). Second, the two state short period 
approximation is valid only for small changes in the velocity from trim. 

Thus, the input for parameter estimation must be designed so that the 
velocity excursions from trim are small in order to ensure validity of a two 
state short period model structure. The velocity is considered an output 
variable with the amplitude constraint (5-18). The model for the velocity (top 
row of (5-16)) must be included in the closed loop model (5-27) used for the 
optimal input design, so that maximum velocity changes from the trim value 
can be limited. Then the two state short period model can be used without 
significant degradation in parameter accuracy due to an invalid model 
structure. In general, using the open loop model (5-16) for inputs which do 
not excite the velocity very much (such as the simple doublet input) produces 
poor estimates of the parameters in the top row of equation (5-16). However, 
the model structure in (5-16) is required to build the a priori closed loop 
model (5-27) needed for the optimal input design. 
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With the velocity perturbations limited (equation (5-18)), all 
subsequent parameter estimation can be based on the two state short period 
model: 
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A longitudinal stick input design proposed by Klein et al [16] was 
applied to the F-18 simulation test aircraft, with the input amplitude adjusted 
by trial and error so that the velocity and angle of attack amplitude constraints 
were satisfied. This input wili be referred to as the compound doublet input. 
The time histories from the experiment using the compound doublet input 
are shown by the solid lines in figure 15. The dotted lines in figures 15(b) and 
15(c) show the model fits for a and q, respectively, from maximum likelihood 
estimation using the open loop model (5-33). Closed loop stabilator deflection 
from the F-18 nonlinear simulation is shown in figure 15(d). 

The optimal input technique was applied using the closed loop model 
(5-27) with the model parameters from Table 6 as the a priori model. The 
optimal input design and time histories from the experiment using the F-18 
simulation test aircraft are shown in figure 16. A fixed time input design was 
produced by setting all Cramer-Rao bound goals to zero, and setting the 
maximum test time to 8.0 seconds. A first order lag with time constant 0.05 
second was used for the longitudinal stick control dynamics. The 0.75 inch 
maximum input amplitude for the optimal input was determined 
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automatically by the algorithm. Maximum input amplitude was initially set 
at 1.0 inch, to agree with the maximum input amplitude for the simple 
doublet. Input amplitude decrements were set at 0.25 inch to design an input 
with an amplitude which could be easily input by the pilot. The algorithm 
thus required one iteration on the problem formulation to arrive at the final 
input amplitude of 0.75 inch. Minimum input pulse width was set at 0.5 
second, and the final time required with zero input amplitude was 0.4 second. 

Table 7 gives the results of maximum likelihood estimation using the 
data from each of the input designs. The parameter values listed in the 
second column were obtained using finite differences with the nonlinear F-18 
simulation. The standard errors for each parameter are lowered using the 
optimal input, except for the small parameter Zg h , which is difficult to 
estimate accurately because of its small magnitude. Figure 16(a) shows the 
optimal input, and figures 16(b) and 16(c) give the F-18 simulation responses 
and the model fits using the open loop model (5-33) and maximum 
likelihood parameter estimates. The closed loop stabilator deflection from 
the F-18 simulation is shown in figure 16(d). Figures 16(e) through 16(h) 
compare the time histories from the optimal input experiment with the 
computed time histories from the optimal input design algorithm using the a 
priori closed loop model (5-27) with parameters from Table 6. The time 
histories from the experiment violate the imposed output amplitude 
constraints because the closed loop model estimated from the doublet input 
does not exactly describe the nonlinear response from the F-18 simulation. 

The a priori closed loop model time histories in figures 16(e) through 16(h) 
satisfy the output amplitude constraints (5-18) and (5-19), as required. 

The example presented here demonstrates that the optimal input 
algorithm can be successfully used for optimal input design problems when 
the estimation of open loop parameters is required, but only closed loop data 
are available. For modern fighter aircraft, this situation has become the norm 
rather than the exception. 
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Chapter VI - Summary and Conclusions 


This dissertation describes a new formulation of the optimal input 
problem for aircraft parameter estimation experiments, the development of a 
new solution algorithm for the optimal input design, and the demonstration 
of the capabilities of the new technique. The examples presented are mainly 
for demonstration purposes. 

The optimal input design for aircraft parameter estimation 
experiments was done in the time domain using the principles of dynamic 
programming. The problem was formulated with the objective of achieving 
specified Cramer-Rao bound goals in a time optimal fashion or minimizing 
the Cramer-Rao bounds for a fixed test time. Optimization in the time 
domain using the dynamic programming approach allowed various practical 
aspects of the input design to be incorporated in a straightforward manner. 
Bellman's principle of optimality was enforced so that the designed inputs 
were globally time optimal, subject to the imposed constraints on the input 
form, the output amplitude constraints, the dynamic and measurement 
models, and the discretization of both time and the constrained output 
variables in the formulation of the dynamic programming problem. 

The first example presented in this work demonstrated that the present 
approach to optimal input design for aircraft parameter estimation improves 
on previous solutions by lowering the Cramer-Rao bounds and/or requiring 
a shorter total test time. 

In the second example, the expanded capability of the optimal input 
design technique was exhibited by application to a fourth order multiple 
input system with restrictive output amplitude constraints. The input design 
solution contained provisions to assist in the practical implementation of the 
designed inputs. 

The third example developed a procedure for the use of the optimal 
input design technique when open loop model parameters are of interest, but 
the system must be tested under closed loop control due to safety 
considerations. Conditions for a real flight test of the F-18 fighter aircraft were 
simulated, and the utility of the optimal input design technique was 
demonstrated. 

The capabilities of the new algorithm can be summarized as: 
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1. Global optimum solution for minimum or fixed flight test time 

2. Multiple input design capability 

3. Closed loop model capability 

4. Practically achievable input design 

5. Direct implementation of output amplitude constraints 

6. Input frequency selection in the time domain 

7. Single pass solution 

8. Well-suited to design studies 

The main contributions of this work are the formulation and practical, 
optimal solution of the input design problem for aircraft parameter 
estimation experiments, along with the demonstrations of the improved 
quality of the optimal inputs and the expanded capabilities of the optimal 
input design technique. 

Based on the findings documented in this work, conclusions and 
recommendations for further study are listed below. 

1. The optimal input design technique should be evaluated with 
regard to pilot acceptability as compared with conventional input 
designs. A ground-based simulator could be used for this study. 

2. Application to a free flight wind tunnel model test, drop model test, 
or full scale flight test would verify the predicted capabilities of the 
optimal input design technique, 

3. Other input design problems which could be addressed using the 
input design procedure outlined here include experiments for data 
compatibility check of the flight instrumentation system, and input 
designs to reduce parameter correlation. The latter should be 
straightforward, since the input design technique currently computes 
the full dispersion matrix recursively, and parameter correlation can be 
computed from the off-diagonal terms of this matrix, normalized by 
the appropriate diagonal elements. 
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4. A study could be done to determine the effect on the optimal input 
design of inaccuracy in the a priori model, i.e., the robustness of the 
optimal input design technique. 

5. Studies concerning various trade-offs among experiment design 
parameters could be done. For instance, what are the effects of different 
measurement noise characteristics or different control system 
dynamics on the optimal input design and the resulting achievable 
parameter accuracies ? 

6. Although the optimal input technique was developed for aircraft 
dynamic models, any other physical system which can be described 
with a linear state space model structure could have been the subject of 
the input design for a parameter estimation experiment. 
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Appendix A -=_ The Information Matrix and Hip 
Cra mer-Rao Lower Bound 


In this appendix, the form of the information matrix used to design 
optimal inputs for parameter estimation is derived. Following this, it is 
shown that the inverse of the information matrix is the theoretical lower 
bound on the parameter covariances. This lower bound is called the 
Cramer-Rao lower bound for the parameter covariances. 

Consider the following state, output, and measurement equations : 

x(t) = f ( x(t), u(t), 0 ) , x(0) = 0 ( , 


y(t) = h ( x(t), u(t), 0 ) 
y m (i) = y(i) + “o(i) 

E(u(i)j = 0 , E{ u(i) u T (j)} = R-5; j 


(A-2) 

(A-3) 

(A-4) 


The sequence of measured output vectors is y m (i), i=l,2,3,...,N, where 
N is the total number of sample times. For a given parameter vector 0 in the 
state and output equations (A-l) and (A-2), respectively, assuming each y m (i) 
is a realization of a Gaussian vector random process, the joint conditional 
probability of realizing the observed sequence of measured output vectors is 
[171: 


P<y m (l W2),..,y m (N) l 0) = 


(2k) 2 IR|'2 


N 


C *P { I-|[ymO)-y<i)] T R' , [y m (i)-y(i)] J 


(A-5) 


The object of parameter estimation is to choose a parameter vector 0 for 0 in 
expressions (A-l) and (A-2) such that the conditional probability in (A-5) is 
maximized. Usually the natural logarithm of equation (A-5) is used for 
computational simplicity. Since probabilities are always positive, and the 
natural logarithm function is a monotonically increasing function for 
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positive arguments, the same 0 vector maximizes both the probability in 
(A-5) and its logarithm. The logarithm of the probability in equation (A-5) 
will be denoted by L( 0 ) and is generally referred to as the log likelihood 
function: 


L( 0 ) = In [p(y m (l),ym(2),...,ym(N) I 0)] (A-6) 

u 0 ) = - y[q in(27t) + in i r i ] + X" ^y m (i) -y (i )] T R 1 [ymG>-yG>] 

i = 1 (A -7) 


The information matrix is defined in terms of the log likelihood 
function [17]: 


f QL(0) auo) l 
30 30* | 


(A-8) 


where E{ } denotes the conditional expectation taken over all possible 
realizations of the measurement vector sequence, y m (i)/ i=l,2,3,...,N and 
conditioned on the parameter vector, 0. Using (A-6), equation (A-8) may also 
be written as: 


M = E 


j d ln[p(Y I 0) ] d ln[p(Y I 0)] \ 

\ 00 00 ' j 


(A-9) 


where Y represents the sequence of measurement vectors, y m (i), i=l,2,3,...,N. 

At this point, a small side excursion is made to find an alternate 
expression for the information matrix, M. Let Y denote a particular 
realization of the measurement vector sequence, y m (i), i=l,2,3,...,N. Then, by 
the definition of conditional probability density functions. 



I 0) dY = 1 


(A-10) 
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The notation here implies that the integral is taken over all possible 
realizations of the measurement vector sequence, y m (i), i=l,2,3,...,N. By 
differentiation. 


3p(Y I 9) (y , e) 5( ln[p(Y I 9 )]) 
50 ’ p 90 


Differentiating (A-10) twice with respect to 0, and using (A-ll) each time. 


5 2 ln[p(Y I 8)] + ^ 91n[p(Y I 8)] 9 ln[p(Y I 0)] 

90 90 T 90 90 1 


p(Y I 0) dY = 0 


Or, using the expectation operator notation, and recalling (A-6), 


E f 9L(0) 9L(8) \ = E ( 9 2 L(8) 

| 50 90 1 j 90 90 1 


Invoking (A-8), 


M = E • 5^1 _ e / 5 2 L(0) ) 


50 98 


90 90 t I 


Equation (A-14) contains the alternate expression for the information matrix, 
M. For the system (A-l) - (A-4), the information matrix may be expressed as 
follows by combining (A-3), (A-4), (A -7), and (A-14), and keeping in mind that 
the measurement vector y m (i) is considered a realization of a random vector 
sequence, and is therefore not a function of 0 : 


M = r'M 

i = A 90 J l 90 


(A-l 5) 
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It will now be shown that the inverse of the information matrix is the 
theoretical minimum value for the covariances of the elements in the 

/v 

estimated parameter vector, 0 . The dependence of L(0) on the input is 
through the state, output, and measurement equations, (A-l) - (A-4); see also 
(A-7). 

Let Y again denote any realized sequence of measurement vectors, 
ym(i)/ i=l,2,3,...,N. Then for an unbiased estimator [1], 

E\0(Y) - 0) = I (e(Y) - 0) p(Y I 0) dY = 0 

•'- 00 (A-l 6) 

Differentiating (A-16) with respect to 0 gives 


■I 


I p(Y I 0) dY + I (e(Y) - e) a P - V J, 9) dY = 0 

30 


(A-l 7) 


The first integral on the left in (A-l 7) is the identity matrix by the definition 
of the conditional probability density function. Substituting into the second 
integral on the left from equation (A-ll), 



(e( Y ) - 0)^?-^ 

30 


10 ) 1 ) 


p(Y I 0) dY = I 


(A-l 8) 


The last equation may also be written as 

E{(0(Y)-0) a(ln[ P (Y r ie)]) } =1 

30 


(A-19) 


The following lemma is now required [1]: 
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Lemma A.l Let x and y be two random p-dimensional vectors. Then, 


E(xxT) * E(xyT|[E(yyT)]-’E(yxT) 


(A-20) 


Proof : 

Let Q be an arbitrary non-random p-by-p matrix. Then, since all 
covariances must be non-negative, 

E( (x - Qy)(x- Qy) T } > 0 (A . 2I) 

Expanding the left side of (A-21), 

E( x x T ) - QE( y xT) - E( x y T ) Q T + QE( y y T j Q T > 0 (A . 22 ) 

E( x x T } > QE( y xT) + E( x yT) Q T - QEj y yT) Q T (A . 23) 

Now choose 

Q=E(xy T )[E(yy r )]- 1 (A . 24) 

so that the two terms on the far right of (A-23) cancel. Substituting (A-24) into 
(A-23) and simplifying. 


E( x x T ) S E| x y T } [E{ y y T )] E[ y x T j 


(A-25) 


This ends the proof of the lemma. 


Returning to (A-19), and invoking lemma A.l (Equation (A-25)) with 


0(Y) - 0 = x 


d ( ln[p(Y I 0)]) 
30 


(A-26) 
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gives 


E |(S(v) - e)Cm - e) T ) , , j E j jata^ijw J ] 0 


Or, using (A-9), 


0(Y)-0 > M' 1 


(A-28) 


Equation (A-28) is called the Cramer-Rao inequality and states that the 
minimum expected values of the parameter covariances for an unbiased 
estimator are given by the inverse of the information matrix. This latter 
matrix is also referred to as the Cramer-Rao lower bound for the parameter 
covariances, or the dispersion matrix, D: 


D = M 1 


(A-29) 


When designing the optimal input for the purpose of parameter estimation, 
only the dispersion matrix is considered, since this separates the merits of an 
input with regard to the parameter estimation accuracy from the merits of 
whatever unbiased estimator may be chosen for the parameter estimation. In 
effect, an input form is evaluated "upstream" of the actual parameter 
estimation calculations by looking only at the theoretical minimum 
parameter covariances possible with that input form, as embodied in the 
elements of the dispersion matrix. The dispersion matrix depends only on 
the sensitivities of the output quantities to changes in the parameters. In 
general, the higher these sensitivities, the more accurately the parameters can 
be estimated, since small changes in parameter values will be manifested in 
large changes in the model output for the system. With large values for the 
sensitivities, the parameter values may be determined very accurately by 
matching measured output with output from the proposed model. The 
diagonal elements of the dispersion matrix are of principal interest, since they 
are the theoretical minimum values of the parameter variances. The square 
root of these diagonal elements are the theoretical lower bounds on the 
parameter standard errors. 
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For the system (A-l) - (A-4), the dispersion matrix may be expressed as 
follows by combining (A-15) and (A-29): 


D = 



R ' 1 


dyO ) 1 

dO j 


(A-30) 


No assumptions of linearity in the system equations (A-l) - (A-4) were 
necessary in this development. 
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A ppendix B - Dynamic Programming 


The algorithm for input optimization described and demonstrated in 
this dissertation employed the principles of dynamic programming. The 
purpose of this appendix is to provide an understanding of the fundamental 
concepts of dynamic programming in the form required for the present work. 
These fundamental concepts are completely contained and described in 
references [10] and [11], along with more specialized information concerning 
dynamic programming and the connection between dynamic programming 
and the calculus of variations. 

Dynamic programming is an efficient means of solving sequential 
multi-stage optimization problems. In figure B-l, the dots represent various 
states of some abstract system. Time is represented by horizontal distance left 
to right, with time zero at the leftmost edge. Time is discretized into stages, so 
that the system states exist only at times which are an integer multiple of one 
stage time. The lines between the dots represent a specific control choice, 
which results in the state trajectory starting from the dot representing the 
initial system state and ending at the dot representing the final system state at 
the next time stage. Each possible control is applied as a constant over one 
stage time, so that each possible state transition (connecting line) corresponds 
to a possible control choice at the initial state. State transitions must always 
move left to right, corresponding to increasing time. All state transitions 
shown in figure B-l are associated with a cost, which is shown as a number 
along the state transition line. The objective is to reach one of the target 
states, represented by the column of dots aligned above the final time stage, by 
starting at the initial state and optimizing the sequence of state transitions 
such that the total cost of the travel from the initial state to one of the possible 
final states is a minimum. Since each state transition corresponds to a control 
choice at the current state, the objective may also be considered to be the 
optimization of a sequence of control choices. The optimal sequence of 
control choices is called the optimal policy. The problem described here is 
often referred to as a discrete multi-stage optimization. 

An exhaustive search for the optimal policy would require the 
investigation of 2 n_1 different paths, where n is the number of discrete time 
stages from the initial condition to one of the target states at the final time 
stage, inclusive. The base 2 in the above expression results from the fact that 
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only two choices are possible at any state for any stage. As the number of 
stages (or the number of control choices) increases, the exhaustive search 
quickly becomes prohibitive in terms of computation time. There is an 
alternative method for the computation of the optimal policy, which is 
known as "the principle of optimality", due to Richard Bellman. In the form 
required here, this principle states : 

An optimal sequence of decisions in a multistage decision process problem 
has the property that whatever the final decision and state preceding the 
terminal one, the prior decisions must constitute an optimal sequence of 
decisions leading from the initial state to that state preceding the terminal 
one. 1 

In figure B-2, the optimal policy was computed by moving left to right 
and using the above principle of optimality. The arrows along the state 
trajectory lines indicate the control and state trajectory which represents the 
best path (lowest total cost) to reach the given state from the immediately 
preceding stage. The italic number at each state indicates the value of the cost 
required to reach that state at that time stage, using the optimal policy 
indicated by the arrows. Figure B-2 represents the results of applying the 
principle of optimality to each state of each stage in a sequential manner, 
moving left to right. The bold line in figure B’2 indicates the overall optimal 
policy and thus the optimal state trajectory, since this path gives the lowest 
cost from the initial state at the initial stage to any of the states at the final 
stage, using an optimal policy. 

The principle of optimality not only gives a systematic method for 
choosing an optimal policy from a large number of possibilities, but also does 
so in a way which realizes a significant savings in computation over an 
exhaustive search. It is this latter quality which makes dynamic 
programming a practical tool for optimization. 


1 Dreyfus, S. E. 1965. Dynamic Programming and the Calculus of Variations . 
New York : Academic Press, p. 8. 
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In Table B-l, the computations required to solve this multi-stage 
optimization problem are compared for the dynamic programming solution 
method versus an exhaustive enumeration, using an increasing number of 
stages. The general expressions for the number of additions and comparisons 
using dynamic programming and exhaustive enumeration are shown at the 
top of Table B-l as a function of the number of stages, n. The total number of 
additions and comparisons required for a dynamic programming solution 
with the grid geometry of figure B-l is on the order of n 2 and 3n, respectively; 
whereas exhaustive enumeration calculations are on the order of n2 and n, 
respectively. Dynamic programming enjoys a considerable computational 
advantage over exhaustive enumeration, even for a relatively small number 
of total stages, as shown in Table B-l. This computational advantage becomes 
even more dramatic when there are more than two control choices at any 
state. 

The grid geometry in figure B-l and B-2 might be considered a discrete 
multi-stage optimization with no state constraints, since for all states at all 
stages except the last, the full range of controls and accompanying state 
transitions are available. That is, a state transition "up" or "down" diagonally 
was possible from any state at any stage except the last. State constraints can 
be introduced by arranging the grid geometry as shown in figure B-3. Here, 
states on the border of some interior stages have only one control choice or 
state transition possibility to the next stage. This can be interpreted as a 
constraint on deviation of the state from the initial state in the abstract sense 
corresponding to vertical displacement within the grid. 

Although the entire discussion here considered the state of some 
abstract system, the output of the same abstract system could just have well 
been used with absolutely no change in any of the development. In fact, the 
present work considers the dots within each grid as representative of some 
abstract system output at the given stage. The reason for this is so that the 
constraints will be on the output variable excursions from their initial values. 

Real dynamical systems are generally modelled as continuous systems 
whose states, outputs, and controls exist at any time. In order to use the 
dynamic programming optimization technique as presented here, it is 
necessary to discretize the continuous problem corresponding to physical 
reality. This is done by dividing an allowable range of output variable 
excursions into small subdivisions. For the case of two constrained output 
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variables, the allowable output space at any stage can be represented by a plane 
region divided into discrete output space boxes, as in figure B-4. Outputs are 
computed from the continuous dynamic and output equations. Those 
outputs which fall inside or on part of the border of a given box in output 
space are considered the same output "state", corresponding directly with the 
dots in figures B-l and B-2. In this way the continuous problem is converted 
into a discrete problem. All controls are constrained to remain constant over 
the time period corresponding to one stage time, so that time can be 
discretized into stages. The cost computed for each control possibility in the 
optimal input design algorithm is directly analogous to the costs associated 
with each state transition in figures B-l and B-2, the latter being represented 
by the numbers next to each state transition line. This completes the 
conversion of the continuous problem into a discrete multi-stage 
optimization, which can be handled using the principle of optimality in the 
manner described above. 
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A ppendix C - Algorithm Flowchart 

The following pages contain a functional flow chart for the algorithm 
used to design optimal inputs for aircraft parameter estimation experiments. 
The numbers in parentheses refer to equation numbers in the text. 
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Figure 1(a) - Square Wave, Energy = 36 
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Figure 1(c) — Square Wave, Energy *= 36 
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Figure 1 (d) — Square Wave, Energy = 1 44 
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Figure 3 - Dynamic Programming Grids 
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Figure 5 (a) - Mehra Input 
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Figure 6 (a) - Chen Input 
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Figure .6 (b) — Chen Input 
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Figure 6 (c) - Chen Input 
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Figure 8 (b) — Optimal Input, Mehra Amplitude 
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Figure 8 (c) — Optimal Input, Mehra Amplitude 
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Figure 10 (c) - Doublet Inputs 
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Figure 10 (d) - Doublet Inputs 
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Figure 10(f) - Doublet Inputs 
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Figure 1 1 (b) - Optimal Inputs 
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Figure 1 1 (c) - Optimal Inputs 
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(d) — Optimal Inputs 
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Figure 1 1 (f) — Optimal Inputs 
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Figure 12 (b) — Optimal Inputs, Increased Amplitude 
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Figure 12 (d) — Optimol Inputs, Increased Amplitude 
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Figure 12 (f) — Optimal Inputs, Increased Amplitude 


122 



AILERON AMPLITUDE (rad) 


12 



123 




RUDDER AMPLITUDE (rad) 


.12 



Figure 13 (b) — Optimal Inputs, Increased Amp., No Control Sequencing 
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Figure 13 (d) - Optimal Inputs, Increased Amp., No Control Sequencing 
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Figure 1 4 (b) Doublet Input 
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Figure 14 (c) Doublet Input 
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Figure 14 (d) Doublet Input 
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Figure 14 (e) Doublet Input 
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Figure 1 4 (f) Doublet Input 
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Figure 16 (c) Optimal Input 
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Figure 16 (e) Optimal Input 
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Figure 16 (h) Optimal Input 


149 






10 




Figure B-2 - Discrete Dynamic Programming Solution 
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Figure B-4 - Continuous to Discrete Grid Conversion 
( Figure 3, repeated ) 
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Xafe lfi l - Saved Variable List 


SflY£.d^ariabIes ,i :u gentTime Stage Only. All Reachable Output Space Boxes 


1.) state vector, x 

(n x 1 real vector) 

2.) state sensitivities, dx/dO 

(n x p real matrix) 

3.) input amplitudes, UASTR 

(m x 1 real array) 

4.) input information, IREPU (3x1 integer array) 

contains: 

a. ) index of last control (integer) 

b. ) number of stages a.) was applied (integer) 

c. ) index of last control different from a.) (integer) 

5.) [d jk ] V j>k. 

(p(p+l))/2 real elements) 

where Dj = [djk] j, 

j=l,2,...,p 

k=l,2,...,p 



i=time stage index 


Saved Variables, Initial i Q - Cu r rent T ime Stage, All Reachable Output Space Boxes 


1.) input sequence 
control indices , IOPTU 


(i x 1 integer vector) 


Table. 2 - Sin gle Input Design Results 


Parameter 

Parameter 

Values 

Cramer-Rao bounds 

Figure 5 
Mehra 
(Ref. 5) 

Figure 6 
Chen 
(Ref. 9) 

Figure 7 
Optimal 
Input, 
Chen 
Amp. 

Figure 8 
Optimal 
Input, 
Mehra 
Amp. 

Figure 9 
Optimal 
Input, 
Mehra 
Amp. 
with lae 

Za 

mm 

HR 

0.0332 

IRB 

Hi 

RRRi 

z «. 

0.005 

Bm 

0.0212 

■ 

0.0214 

0.0213 

Ma 


HUB 


RHB 

■ 

0.0647 



Hi 


0.1167 

m 

0.1053 

M s e 

■i 

HR! 

0.0740 

■ 


0.0559 

Max.Amp. 

(rad.) 

12.5 

8.792 

8.792 

12.5 

12.5 

Total Time 
(sec.) 

4.0 

3.96 

— 

3.68 

3.04 

3.20 
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labk 3 - Single Input Maximum Likelihood Results 



True 

Value 

Mehra 
(Ref. 5) 

Chen 
(Ref.9 ) 

Optimal 

Input, 

Chen 

Amp. 

Optimal 

Input, 

Mehra 

Amp. 

Optimal 

Input, 

Mehra Amp. 
with lag 


B5S8S 

mam 

n 


BKII99 

WBWSWW 

E9 

WSBL 

■nm 

0.0329 

BSK11 

0.0362 

0.0344 

Z «e 

m 

■ 

0.0227 

0.0232 

Bl 

0.0238 



0.0261 

0.0209 

0.0232 


wmm 

M a 

m 

Bum 

BIB 

■I 

IBH 

warn 



BUSH 

0.0651 

mm 


mrnm 

B 

-1.5880 

WM 

_ -1..5921 

-1.5979 


IHMBI 

ESMi 


0.1829 

wm 

0.1158 

0.0978 


M 6e 


HU 

■ 

m 

H 

wm 

°M 8 , 

■■ 

HfflHI 

RB9 

■m 

BillSi 
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Table 5 - F-18 Geometry and Mass Characteristics 


Total length, m 

17.07 

Wing: 

Are?, m 2 

37.16 

Span, m 

11.41 

Mean geometric chord, m 

3.51 

Aspect Ratio 

3.5 

Quarter-chord sweep angle, deg 

20.0 

Horizontal Tail: 

Area (wetted), m 2 

16.35 

Span, m 

6.58 

Mean geometric chord, m 

1.91 

Aspect ratio 

2.4 

Quarter-chord sweep angle, deg 

42.8 

Moment arm (c.g. at 0.25 m.a.c.), m 

5.12 

Weight, lbs 

31,748.9 

Inertia: 

Ixx/ slug-ft 2 

22,294.7 

lyy, slug-ft 2 

123,095.3 

I ZZ/ slug-ft 2 

138,116.7 

Ixz, slug-ft 2 

-1797.4 

c.g. location: 

fuselage station, in 

456.0 

buttock line, in 

0.0 

water line, in 

103.4 
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Table B-l - Computational Load Contrast 


Dynamic Programming Exhaustive Enumeration 



Additions 


Additions 

Comparisons 

n 

n-1 

£ 2i = n 2 -n-2 
1-2 

2(n-2) + (n-l) = 3n-5 

2"' 1 * (n-2) 

n-1 

H 

18 

10 

48 

4 

10 

70 

25 

4096 

9 

H 

340 

55 


18 
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